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ABSTRACT 

N-body simulations have demonstrated a correlation between the properties of 
haloes and their environment. In this paper, we assess whether the ellipsoidal collapse 
model, whose dynamics includes the tidal shear, can produce a similar dependence. 
First, we explore the statistical correlation that originates from Gaussian initial condi- 
tions. We derive analytic expressions for a number of joint statistics of the shear tensor 
and estimate the sensitivity of the local characteristics of the shear to the global ge- 
ometry of the large scale environment. Next, we concentrate on the dynamical aspect 
of the environmental dependence using a simplified model that takes into account the 
interaction between a collapsing halo and its environment. We find that the tidal force 
exerted by the surrounding mass distribution alters the axes collapse and causes haloes 
embedded in overdense regions to virialize earlier. The environment density is the key 
parameter in determining the virialization redshift, while the environment asphericity 
primarily contributes to the increase in the scatter of the critical collapse density. An 
effective density threshold whose shape depends on the large scale density provides a 
good description of this environmental effect. Such an interpretation has the advan- 
tage that the excursion set formalism can be applied to quantify the environmental 
dependence of halo properties. We show that, using this approach, a correlation be- 
tween formation redshift, large scale bias and environment density naturally arises. 
The strength of the effect is comparable, albeit smaller, to that seen in simulations. 
It is largest for low mass haloes (M <c M*), and decreases as one goes to higher 
mass objects (M > M*). Furthermore, haloes that formed early are substantially 
more clustered than those that assembled recently. On the other hand, our analytic 
model predicts a decrease in median formation redshift with increasing environment 
density, in disagreement with the trend detected in overdense regions. However, our 
results appear consistent with the behaviour inferred in relatively underdense regions. 
We argue that the ellipsoidal collapse model may apply in low density environments 
where nonlinear effects are negligible. 
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1 INTRODUCTION 

In standard scenarios of structure formation, dark matter 
haloes grow hierarchically from initially small, Gaussian 
fluctuations. Properties of haloes can be studied in great de- 
tail using both N-body simulations and analytic models. The 
remarkably useful Extended Press- Schechter (EPS) theory 
predicts halo mass functions (Press & Schechter 1974; Bond 
et al. 1991), merging histories (Lacey & Cole 1993; Sheth 
& Lemson 1999b; Van den Bosch 2002; Neistein, Van den 
Bosch & Dekel 2006) and spatial clustering (Mo & White 
1996; Mo, Jing & White 1997; Catelan et al. 1998; Sheth & 
Lemson 1999a) that are in reasonable agreement with the 
simulations. This analytic approach is based on the spheri- 
cal collapse model (Gunn & Gott 1972) . In this Lagrangian 



approximation, haloes are identified in the initial conditions 
and a single parameter, the initial density contrast, is needed 
to characterise their epoch of formation (Press & Schechter 
1974). The collapse of a halo occurs when the linear density 
reaches a critical threshold. The fundamental properties of 
dark matter haloes are then obtained from the statistics of 
trajectories of the linear density field as a function of the 
smoothing scale (e.g. Bower 1991; Bond et al. 1991; Lacey 
& Cole 1993; Kauffmann & White 1993; Kitayama & Suto 
1996; Sheth & van de Weygaert 2004). Although this spheri- 
cal approximation works well until the first orbit crossing, it 
may not be accurate since perturbations in Gaussian density 
fields are inherently triaxial (Doroshkevich 1970; Bardeen et 
al. 1986; Jing & Suto 2002). Furthermore, the initial shear 



© 0000 RAS 



2 Vincent Desjacques 



field rather than the density has been shown to play a cru- 
cial role in the formation of nonlinear structures (e.g. Hoff- 
man 1986, 1988; Peebles 1990; Dubinski 1992; Bertschinger 
& Jain 1994; Audit & Alimi 1996; Audit, Teyssier & Alimi 
1997). 

Unlike the spherical model whose dynamics depends on 
a single parameter only (the density), the ellipsoidal model 
that follows the evolution of triaxial perturbations can be 
used to ascertain the influence of the (external) tidal shear 
on the properties of collapsed regions. The gravitational col- 
lapse of homogeneous ellipsoids has been investigated by 
numerous authors over the past decades (e.g. Lynden-Bell 
1964; Lin, Mestel & Shu 1965; Fujimoto 1968; Zeldovich 
1970; Icke 1973; White & Silk 1979; Barrow & Silk 1981; 
Lemson 1993; Eisenstein & Loeb 1995; Hui & Bertschinger 
1996). In the formulation of Bond & Myers (1996), initial 
conditions and external tides are chosen to recover the Zel- 
dovich approximation in the linear regime. The dynamic of 
ellipsoidal collapse can be incorporated in various ways in 
the Press-Schechter formalism to predict the properties of 
haloes (Monaco 1995, 1997a, 1997b; Lee & Shandarin 1998; 
Chiueh & Lee 2001; Sheth & Tormen 2002). As pointed 
out by Sheth, Mo & Tormen (2001), the inclusion of non- 
sphericity in the dynamics introduces a simple dependence 
of the critical collapse density on the halo mass. The result- 
ing first crossing distribution yields a better fit to the halo 
mass functions measured in N-body simulation (Sheth & 
Tormen 2002). However, other modifications to the original 
excursion set approach, such as the inclusion of non-radial 
degrees of freedom, might also improve the theoretical mass 
function (Audit et al. 1997; Del Popolo & Gambera 1998). 
It would thus be very desirable to identify additional dis- 
tinctive predictions of the ellipsoidal collapse model beyond 
the mass function and bias to further test this theory. 

While earlier numerical studies have not provided any 
conclusive evidence for a dependence of halo properties on 
environment (Lemson & Kauffmann 1999; Percival et al. 
2003; Zentner et al. 2005), recent numerical investigations 
indicate that, at fixed halo mass, haloes in dense regions 
form at (slightly) higher redshift than in low density environ- 
ments (Sheth & Tormen 2004; Avila-Reese 2005; Harker et 
al. 2005). Using the Millennium Run (Springel et al. 2005), 
Gao, Springel & White (2005) have convincingly shown that 
the clustering of haloes of a fixed mass depends on forma- 
tion time. This dependence is strong for haloes with mass 
less than the typical collapsing mass A/*, and fades rapidly 
for M > Af*. Subsequent studies have demonstrated that 
many other halo properties, such as spin parameter or con- 
centration, correlate with the halo assembly history (Maul- 
betsch et al. 2006; Wechsler et al. 2006; Zhu et al. 2006; Gao 
& White 2007; Wetzel et al. 2007; Jing, Suto & Mo 2007). 
Also, haloes that have undergone major mergers may be 
more strongly clustered relative to other haloes of the same 
mass (e.g. Furlanetto & Kamionkowski 2006). These results 
call into question the simplest descriptions of structure for- 
mation based on the statistics of random walks (Bond et al. 
1991; Lacey & Cole 1993; White 1996). However, relaxing 
the assumption of sphericity and/or sharp fc-space filtering 
can introduce a dependence on environment (Bond et al. 
1991; White 1996; Sandvik et al. 2007). In the ellipsoidal 
model, the time required for a given overdensity to virialize 
increases monotonically with the initial shear (Sheth, Mo 



& Tormen 2001). Therefore, as recognised by Wang, Mo & 
Jing (2006), the ellipsoidal dynamics should give an environ- 
mental effect owing to the tidal field generated by the large 
scale environment. 

In the present paper, we take an analytic approach and 
assess whether the ellipsoidal collapse model can produce 
an environmental dependence (also termed 'assembly bias') 
similar to that seen in N-body simulations. We start with the 
statistical dependence that arises in correlated (Gaussian) 
initial conditions. We extend the results of Doroshkevich 
(1970) to the joint statistics of the shear tensor. We derive 
conditional distributions and quantify the extent to which 
the asymmetry of initially triaxial perturbations is sensitive 
to the geometry of the large scale environment. Next, we 
investigate the dynamical aspect of the environmental de- 
pendence using a simplified model that takes into account 
the interaction between a triaxial protohalo and its environ- 
ment. We find that the tidal force exerted by the surround- 
ing mass distribution affects the axes collapse and causes 
haloes embedded in large overdensities to virialize earlier. 
A moving barrier whose shape depends on the environment 
density provides a good description of this environmental ef- 
fect. This enables us to apply the EPS formalism in order to 
estimate the environmental dependence of halo properties. 
Our approach thus is very different than the multidimen- 
sional extension presented in Sandvik et al. (2007). 

The paper is organised as follows. Section [5] briefly re- 
views the basic concepts associated with the ellipsoidal col- 
lapse model. Section Jj3]is devoted to the statistical correla- 
tion between the local properties of the shear and the large 
scale environment (Appendix [B] details a delicate step of 
the calculation). |4] investigates the dynamical origin of the 
environmental dependence, focusing on the distribution of 
halo formation redshift, large scale bias and alignment of 
spin parameter. Non-Markovianess and tidal interactions as 
potential sources of environmental dependence are discussed 
in Sj5] A final section summarises our results. 



2 THEORETICAL BACKGROUND 

We emphasise the role played by the shear in current theories 
of structure formation and introduce the basic definitions 
and relations relevant to the study of environmental effects 
in the ellipsoidal collapse model. 

2.1 Shear tensor 

The comoving Eulerian position of a particle can be gen- 
erally expressed as a mapping x — q + S(q,t), where 
q is the Lagrangian (initial) position and S is the dis- 
placement field. In the Zeldovich approximation (1970), 
the displacement field is S(q,t) = —D(t)X7<&(q), where 
= <t>{q,t) / A-KGpm{t)a 2 D{t) is the perturbation poten- 
tial (<f)(q,t) is the Newtonian gravitational potential), p m is 
the average matter density and D(t) is the linear growth 
factor (Peebles 1980). The second derivatives of the per- 
turbation potential define the deformation tensor (or strain 
field) Dij = didjQ. For convenience, we introduce the real, 
symmetric tensor 
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where a = a(R) is the rms variance of density fluctuations 
smoothed on scale R (see EI3.1I below). We will henceforth 
refer to £jj as the shear tensor. Let Ai > A2 > A3 designate 
the ordered eigenvalues of £y. An important quantity is the 
probability distribution of the ordered set (Ai,A2,As), first 
derived for Gaussian random fields by Doroshkevich (1970), 



P(Ai,A 2 ,A 3 ) 



where 



15 J 
8-kVE 



A(A)e 



-8s?(A) + *f s 2 (A) 



A( a; )=det(xf- J ) = J] (^-^) 



(2) 



(3) 



l<i<i<iV 



is the Vandermonde determinant in the arguments x% , i = 
1, . . . , N, and s n (x) are the elementary symmetric functions 
of degree n (Weyl 1948). For three variables Xi,x%,X3, the 
first few elementary functions are 



Si(x) = Xl + X-2 + X 3 

S2(x) = X1X2 + X1X3 + X2X3 

S 3 (x) = X1X2X3 . 



(4) 



If xi, . . . ,xn are the eigenvalues of a matrix X, the func- 
tions s„(x) can be written in terms of the traces of power 
of X, trX fc with k, I = 0,1,.... For instance, S2(x) = 



(1/2) [(trX) 2 - tr (X 2 )] etc. 



2.2 Geometry of the initial density field 

As shown by Bond, Kofman & Pogosyan (1996), the fila- 
mentary pattern seen in N-body simulations (e.g. Park 1990; 
Bertschinger & Gelb 1991; Cen & Ostriker 1993; Springel et 
al. 2005 for a recent example) is a consequence of the initial 
spatial coherence of the shear tensor. In this Cosmic Web 
paradigm, the correspondence between large scale structures 
in the evolved density field and local properties of the shear 
tensor in the initial conditions, and the knowledge of the 
probability P(Ai,A2,As), allows us estimate the morphol- 
ogy of the large scale matter distribution. 

The geometry of the primeval density field depends on 
the signature of the ordered sequence of shear eigenvalues. 
If, in a given region, the largest eigenvalue only is positive 
(H ), there is contraction along one direction and ex- 
pansion in the other two so that a pancake will form. If 
two eigenvalues are positive while the third one is nega- 
tive (+ H ), collapse occurs along two directions and a 

filament will form. The probability for these two configura- 
tions, ~ 0.84, is much larger than the probability that all 
three eigenvalues are positive, P(+ + +) — 0.08. However, 
these values depend strongly on the density enhancement 
of the region under consideration. While filaments or sheet- 
like configurations are favoured when v <l-5, one encoun- 
ters predominantly spherical-like mass concentrations above 
v ~ 1.5 (e.g. Bernardeau 1994; Pogosyan et al. 1998). Note 
that several methods have recently been proposed to quan- 
tify precisely the geometry and topology of density fields 
(e.g. Sahni et al. 1998; Colombi, Pogosyan & Souradeep 
2000; Hanami 2001; Novikov, Colombi & Dore 2006; Gleser 
et al. 2006; Aragon-Calvo et al. 2007a). 

The sequence in which these structures form is still an 
open problem. We briefly note that, in the pancake pic- 
ture, the collapse proceeds in the order pancakes-filaments- 



-clusters (Lin, Mestel & Shu 1965; Zeldovich 1970; Arnold, 
Shandarin & Zeldovich 1982). Supporting evidence comes 
from several N-body simulations showing first pancake-like 
collapse (e.g. Shandarin et al. 1995). 

Identifying the precursors of haloes (protohaloes) in the 
initial conditions is another unresolved issue, despite the 
major advance made in the analysis of (Gaussian) random 
fields (Doroshkevich 1970; Adler 1981; Peacock & Heavens 
1985; Bardeen et al. 1986; Bond & Myers 1996). Shandarin & 
Klypin (1984) have shown by means of simulations that mas- 
sive clusters with M ;>10 15 Mg/A are initially located close 
to local maxima of the smallest eigenvalue, with A3 > 0. On 
the other hand, recent N-body investigation (Porciani, Dekel 
& Hoffman 2002a,b) indicate that a large fraction of small 
mass haloes M < M+ are rather associated with primeval 
configurations of signature (+H — ). Since, at present, there 
is no reliable alternative to the Press-Schechter prescription, 
we assume that haloes form out of initially spherical patches 
of size R, so that the relevant averages are taken over spher- 
ical Lagrangian regions. This assumption is quite unrealistic 
in the light of numerical results, but it greatly facilitates the 
calculation. We shall also adopt the usual critical density cri- 
terion issued from the spherical collapse (see below). This 
restrictive collapse condition does not guarantee a strictly 
positive signature (+ + +)• But in the ellipsoidal collapse 
model, the formation of a bound object can occur even if 
A3 < : once the shortest axis has collapsed, the nonlinear 
density causes the other two axes to collapse very rapidly 
(Bond & Myers 1996). 

2.3 Ellipticity, prolateness and critical collapse 
density 

The eigenvalues Ai can be equivalently parametrised in 
terms of the shear ellipticity e and prolateness p, where 

Ai — A3 Ai — 2A2 + A3 



2;/ 



2v 



(•») 



and v — S/a = Ai + A2 + A3 is the density contrast of the 
region under consideration. The ordering constraint implies 
that e > if v > and e < if v < 0. In all case, the 
shear prolateness is — e < p < e. In this parametrisation, 
extreme sheet-like (oblate) structures have p <; + e while 
extreme filaments (prolate) have p ;> — e. Doroschkevich's 
formula can be used to write down the distribution g{e,p\u) 
of ellipticity e and prolateness p for a given density contrast 
v (e.g. Bardeen et al. 1986), 



, 1125 i 2 

g{e,p\v) = —=== e (e 



i, 2 (3 e 2 +p 2 ) 



(6) 



For all v, the maximum of this distribution occurs at 

e m (y) = l/(VEu), p m (u) = 0, (7) 

while the variances of e and p are (for v > only) 

{e 2 \v)-{e\v) 2 = (197r-54)/(607r^ 2 ), (p 2 \u) = 19/(20^ 2 ) .(8) 

Notice that the most probable value e m is comparable to 
the mean ellipticity (e\v) = 3/(V 10nv) ~ 1.20e m . It de- 
pends on the smoothing scale R through v oc S/a(R). At 
fixed R, denser regions are more likely to be spherical than 
less dense regions while, at fixed v, larger regions are more 
likely to be spherical than smaller ones. The scatter in the 
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asymmetry parameters increases strongly with decreasing 
density and/or halo mass. 

In the excursion set formalism (Bond et al. 1991; Lacey 
& Cole 1993) , the critical collapse density encodes the details 
of the collapse dynamics. In the spherical collapse model, 
the dynamics in a given cosmological background is gov- 
erned by a single parameter, namely the density. A top-hat 
perturbation of overdensity v = 5 sc /a (linearly extrapolated 
to present epoch) collapses at redshift z = 0. The (linear) 
critical density threshold 8 SC depends on the cosmology. We 
have <5 SC = 1.673 for the cosmological parameters considered 
here (Eke, Cole & Frenk 1996; Navarro, Frenk & White 1997; 
Lokas & Hoffman 2001). On the other hand, in the ellipsoidal 
collapse dynamics, the evolution of a perturbation depends 
on the values of e, p and v. The critical density threshold 
8 CC is always larger than the spherical value 8 SC and is very 
sensitive to the initial shear (Sheth, Mo & Tormen 2001). 



3 ENVIRONMENTAL EFFECT FROM THE 
STATISTICS OF THE INITIAL SHEAR 

The statistics of the shear tensor for Gaussian random fields 
has been pioneered by Doroshkevich (1970) to study the 
formation of large scale structures. In his seminal paper, 
Doroshkevich calculated the probability distribution of the 
shear eigenvalues and ascertained the amount of material be- 
ing incorporated in a pancake. Later, Doroshkevich & Shan- 
darin (1978) reexamined the formation of sheet-like struc- 
tures and derived a distribution function for the largest 
eigenvalue of the shear tensor. Recently, Lee & Shandarin 
(1998) computed conditional probability distributions for 
individual shear eigenvalues to obtain an analytic approx- 
imation to the halo mass function. 

Here, we derive expressions for a number of joint statis- 
tics of the shear tensor upon the assumption of Gaussianity. 
This enables us to quantify the importance of the statisti- 
cal correlation between the asphericity of triaxial collapsing 
regions and the shape of their large scale environment. 

3.1 Analytic considerations 

We confine our calculation to the case in which the com- 
ponents ^ij(x) and £,ki(x') are smoothed on different scales, 
but the joint distribution is evaluated at a single comoving 
position x = x' . This is most relevant to the issues consid- 
ered in this paper. The central results of this Section are the 
joint probability distribution of shear eigenvalues, eq. (I23p . 
and the conditional probability for the shear ellipticity and 
prolateness, eq. (|27[) . 

3.1.1 Spectral parameter 

We begin with the two-point correlation functions of the 
shear tensor. The general form of these correlations is given 
in Appendix SjX] Evaluated at a single comoving position x, 
they take the simple form 

{£,ij(x)£,ki(x)} = — (SijSki + StkSji + Su&jk) , (9) 

10 

where £y and £kl are smoothed on comoving scale Rq and 
Ri respectively. The spectral parameter 
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Figure 1. The spectral parameter 7 as a function of the ra- 
tio Rq/Ri for three different smoothing lengths Ri = 0.5, 2 
and 6 h~ 1 Mpc (solid curves from top to bottom). These val- 
ues correspond to a mass scale Mi = 3.3 X 10 10 , 2.2 X 10 12 and 
5.8 X 10 13 Mq//i, respectively. The dotted curves show the scaling 
7 oc (Ri/Ro)~ ( 3 + n eff)/ 2 ( see text). In all cases, 7 is large even 
on scales Ro 2> Ri ■ 

1 f°° 

-/= / dlnfcA^(fc) W(Ra,k)W(Ri,k) , (10) 

ffoffi J 

< 7 < 1, is a measure of the correlation between these 
scales. Here, Af(fc) = k 3 P s (k)/2n 2 is the dimensionless, 
linear density power spectrum (Peebles 1980) and ao and 
a\ are the rms variances of density fluctuations smoothed 
on scale Ro and Ri, respectively. Ri is the comoving char- 
acteristic scale of the spherically symmetric window func- 
tion W(Ri,k). Many choices are possible for this filtering 
function. We will adopt a top-hat filter throughout this pa- 
per. The top-hat smoothing radius Ri defines a mass scale 
Mi — (47r/3) /9 m i? 3 so that, for a given power spectrum, 
(Ti, Mi and -Ri are equivalent variables. Note also that, in- 
stead of writing explicitely the smoothing radius, we will 
use subscripts to distinguish quantities at different smooth- 
ing lengths. We will reserve the subscript 1 for haloes and 
the subscript for the environment, assumed uncollapsed at 

2 = 0. 

Fig. [T] shows the correlation strength 7 as a function of 
Ro for three different values of Ri : 0.5, 2 and 6 /i -1 Mpc 
(curves from top to bottom). The curves have been com- 
puted for a ACDM model of spectral index n s — 0.96 
and normalisation as = 0.83 using the fitting formulae 
of Eisenstein & Hu (1999). This choice is consistent with 
the constraints inferred from the latest CMB measurements 
(WMAP3, see Spergel et al. 2006). For the special case of 
a power-law spectrum Ps(k) oc k", the parameter 7 scales 
with the smoothing lengths Ro and R\ as 

/ Ri \ < 3 +")/ 2 

Hf) ■ < u > 

Strictly speaking, this expression is valid for a power-law 
spectrum only. However, it provides a reasonable approxi- 
mation in the ACDM cosmology considered here if the spec- 
tral index n is replaced by an effective index n a ff(k) = 
dlnP^(fc)/dlnfc evaluated on scale k ~ 1/Ro- The scaling 
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7 oc (R 1 /R )- {3+ "'= a)/2 is plotted as dotted curves in Fig.[U 
Recall that the spectral index is close to n e « ~ -2 on co- 
moving scales 1 — 10 ft - Mpc. 

3.1.2 Joint statistics of the shear tensor 

Owing to the symmetry of £jj, only six components are in- 
dependent. We adopt the notation of Bardeen et al. (1986) 
and label them by £a, where the components A = 1, . . . , 6 
of the six-dimensional vector refer to the components ij = 
11, 22, 33, 12, 13, 23 of the tensor. The joint probability dis- 
tribution P(£o,£i) of the shear tensor £o and £i, smoothed 
on scale Ro and Ri respectively, is given by a multivari- 
ate Gaussian whose covariance matrix M has 12 dimen- 
sions. This 12 x 12 matrix may be partitioned into four 
6x6 block matrices, Mi = (£i£i ) in the top left corner, 
M2 = (£o£o } m the bottom right corner, B = {£, £,J} m the 
bottom left corner and its transpose B T in the top right cor- 
ner. Following Bardeen et al. (1986), we transform the six 
dimensions {£o,A! £i,A, A = 1, 2, 3} to a new set of variables 
{u k ,v k ,w k ,k = 0, 1}, where 

Uk = £k,l + £fc,2 + £fc,3 
1 



Vk 



(£fc,i — £it,3) 



ui k = - (£*.,! - 2 f; fc , 2 + £4,3) 
With these definitions. 



(12) 



= 1, («k) = — , (w 2 k ) = - 

(u wi) = 7; (vovi) = — 1 {w wi) = - , (13) 

The other correlations are zero. For the six remaining com- 
ponents {£o,a, £i,A, A = 4, 5, 6}, the correlations functions 
are 

1 



(£o,a£o,b) = (£o,a£o,b) 



15 



5ai 



(£o,a£i,b) 



1_ 

15 



(14) 



All the cross-correlations between u k , v k , w k and £o,a, £i,a 
vanish. The block matrices Mi, M2 and B are diagonal in 
the basis introduced above, 



Mi = M 2 = 



C 
1/15 



B = 



7 C 

71/15 



(15) 



where C = diag(l, 1/15, 1/5) and I is the 3x3 identity ma- 
trix. The quadratic form which appears in the joint proba- 
bility distribution 

1 



P(£o,£i)d£ d£i 



(2 7 r) 6 |detM| 1 /2 



e~ Q(&,Sl) d£ d£i 



(16) 



where detM is the determinant of the covariance matrix M, 
can be computed easily using Schur's identities. The result 
may be expressed in terms of the elementary symmetric 
functions Q or, equivalently, in terms of traces, 

Q(£o,£i) = ! 



iy{5[tr (£ 2 )+tr (£?) -2 7 tr(£ £i)] 



4(1-7*) 
[(tr£ ) 2 + (tr£i) 2 - 2 7 (tr£ )(tr£i)]} 



(17) 



and the square-root of the determinant is given by 
|detM| 1/2 = (20/15 6 )(l - 7 2 ) 3 . The results are described 



by the correlation strength 7 solely. The invariance under 
rotation P(£o,£i) = P (R£oR T , R£iR T ) , where R is a real 
orthogonal symmetric 3x3 matrix, requires that P be a 
symmetric function of the eigenvalues, and thus a function 
of tr (£o£i) > k, I = 0, 1, . . ., regardless the statistical proper- 
ties of £ij . The expression (|17p follows from our assumption 
of Gaussianity. Note also that no assumptions have been 
made so far about the coordinates. 



3.1.3 Joint distribution of the eigenvalues 

Lee & Shandarin (1998, see their Appendix B) have com- 
puted the joint probability distribution of the eigenvalues of 
the deformation tensor for a sharp fc-space filter. However, 
they assume that both principal axis frames are aligned, 
which is not true in general. 

To obtain the joint probability distribution of the or- 
dered eigenvalues of the shear tensor, we choose a coordinate 
system such that the coordinate axes are aligned with the 
principal axes of £o- Let a and A be the diagonal matrices 
consisting of the three ordered eigenvalues ai > 02 > Q3 
and Ai > A 2 > A3 of the deformation tensors £0 and £1, 
respectively. The principal axis are labelled according to 
this ordering. With this choice of coordinate, £0 = ct and 
£1 = RAR T , where R is an orthogonal matrix that defines 
the orientation of the eigenvectors of £1 relative to those of 
£0- To preserve the orientation of the principal axis frames, 
we further impose the condition that the determinant of 
R must be +1. Namely, R belongs to the special orthog- 
onal group SO (3). The properties of the trace imply that 
tr£i = trA and tr (£ 2 ) = tr(A 2 ). We note, however, that 

the term tr(£o£i) = tr (RAR T a) depends on the rotation 
matrix. 

The joint probability distribution P(a, A) is obtained by 
integrating over the rotations that define the orientations 
of the orthonormal eigenvectors of £0 and £1. The volume 
measure d£ for the space of real 3x3 symmetric matrices 
can be expressed in terms of the non-increasing sequence of 
eigenvalues ij as 



d£ = 8tt 2 A(i)d 3 tdR . 



(18) 



Here, dR is the invariant measure on the group SO(3) nor- 
malised to J dR = 1, A(t) is the Vandermonde determi- 
nant (eq. [3]l and d 3 t — dtidt2dt3. Since the quadratic form 
Q depends only on the relative orientation of the two or- 
thonormal triads, we can immediately integrate over one of 
the SO(3) manifolds. The relevant volume is 87r 2 /4 = 27r 2 . 
The factor 4 comes from not caring whether the rotated axis 
points in the positive or negative direction (see e.g Bardeen 
et al. 1986). The essential problem is the calculation of the 
integral over the rotations that define relative, distinct triad 
orientations. Appendix !jB] shows that the integral over the 
second SO(3) manifold can be cast into the form 



dR exp [j0tr (RAR T a)] 



W(/3e- 



(19) 



SO(3) 



where the function W depends on 7 through the parameter 
,3(7) = (15/2) 7 /(l — 7 2 ). The four independent variables e+, 
e_, e a and e\ are combinations of the six eigenvalues of a 
and A (Wei & Eichinger 1990). We choose e+ = (1/3) tratrA, 
e_ = (1/3) tratrA, e a — (ai — a 2 )/tra, e\ — (Ai — Aa)/trA, 
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where tra — tra — 3Q3 and trA = trA — 3A3. With this 
parametrisation, —00 < e+ < 00, e_ > and < e a , e\ < 1. 
The function W(/3e_, e a , e\) can be written down as a dou- 
ble integral (eq. |B9j. We have found that a fifth-order ex- 
pansion about /3e„ = (eq. IB12[) is accurate to within 2 per 
cent in the range < /3e_ <l-5. We use this truncated se- 
ries in the computation of the group integral (|19fl . The joint 
probability distribution P(a, A)d 3 a d 3 A may now be formu- 
lated as 



P(a, A)d 3 a d 3 A 



15° 



1 



7 2 )"V(/3 e . 



320tt 2 

xe~ Qol+f3c+ A(q)A(A) d 3 ad 3 A , (20) 



where the quadratic form Qd is a function of a and A, 
3 



- [(tra) 2 + (trA) 2 - 2 7 (tra) (trA)] } 



(21) 



Note that, in the limit 7 <g 1, /3 -> and the joint prob- 
ability distribution P(a, A) tends, as it should be, toward 
the product of the individual one-point probability distribu- 
tion eq. ([21 Using Bayes' theorem, we can easily derive the 
distribution of eigenvalues Xi given the values of at. This 
conditional probability may be written as 



P(A|a) 



15 d 



x exp 



(l- 7 2 ) 3 W(f3e- y e a ,e x )e^A(\) (22) 

trA) 2 " 



157 2 tr (a 2 ) + 15tr (A 2 ) - 3 ( 7 t 



4(l-7 2 ) 



It can be verified by direct numerical integration that the 
probability distribution (|23[) is normalised to unity (we have 
used the multidimensional integrator DCUHRE described in 
Berntsen, Espelid & Genz 1991). 



3.1.4 Joint distribution of the shear ellipticity and 
prolateness 



The results of Section i]3.1.3l can be conveniently expressed 
in terms of the density contrast v, shear ellipticity e and 
shear prolateness p. The joint probability distribution for 
these variables, g (ei,pi, V\, eo,po, ^o)) is readily obtained 
from the following coordinate transformation, 



Qi = y (1 + 3e +po) 



Ai = y (l + 3ei+pi) 



U ° It 
«3 = T (1 



2po) 

3e +Po) 



A 3 = f-d 



2pi) 

3ei +pi) . (23) 



We have, for example, tr(a) = uq, d 3 a = (2/3)^ 2 d^odeodpo 
and A(a) = 2i/q eo (eg — p 2 ,) . With this parametrisation, the 
quadratic form Q01 becomes 



2(l-7 2 ) 



\yl (3eo + pi) + v\ (3ef + pi) + jvqv-l] 



1 

+ 2 



(y x - 7 i/ ) 2 

— 1 2 ^0 

1 — 7-° 



(24) 



The other variables are simply e+ = (1/3)^0^1, £- = 
(l/3)i/o^i (3e -po) (3ei -pi), e a = (e + po) / (3e - po) 
and e A = (ei + pi) / (3ei - pi). 

Let us introduce the notational shorthands = 
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Figure 2. The 20 and 68 per cent contours of the conditional 
probability distribution g(ei,pi\vi, 0) for several values of the 
correlation strength, 7 = 0, 0.4, 0.6 and 0.8 (from top to bottom). 
The values of (eojPOi^o) are typical of a (prolate) filament. The 
dashed lines indicate the boundary of the domain \p\ > e. 



(eo,Po,fo) and 1 = (ei,p%,ui). We wish now to calculate 
the conditional probability g{e\,pi\v\,Q) of having an ellip- 
ticity ei and prolateness pi given a density V\ on scale Ri, 
and given the values (ec^po^o) on scale Ro- This probabil- 
ity will be useful to estimate the effect of the large scale 
environment on the primeval distribution of shear ellipticity 
and prolateness. Bayes' theorem implies that 



g(ei,pi\ui,0) 



5(0,1) 



(25) 



<?K0) 

Since (za£o,a) = 7/3 if A = 1,2,3 and the other cross- 
correlations are zero, the calculation of the denominator is 
straightforward. We have g(vi,0) = <?(i/i|^o) <?(0), where 

^2 " 



5(^1 M = -j= (l - 7 2 ) 1/2 exp 



2?r 



7^0 ) 



2 (l- 7 2) 



(26) 



as the density contrast v is independent of the asymmetry 
parameters e and p. With these informations, the conditional 
probability distribution can be expressed as 



ff(ei,p 2 K,0) = 6(0, l)5(ei,pi|i/i) , 
where 

/ 1 \ 1125 5/2 2 \ 
ff(ei,Pi|^i) = rT7r - V\ei [e 1 -p L ) 



(27) 



(28) 



is the distribution without the environmental constraint 
(see, e.g., Sheth, Mo & Tormen 2001), and 

-5/2 . 



6(0,l) = (l- 7 2 ) ' W(/3e. 



(29) 



x exp 



57 



2(l-7 2 ) 



[vl (3eo + Po) + v \ (3e 2 + p 2 )] 



is a correction factor which results from the constraint "0" . 
As expected, in the limit where the correlation becomes 
weak, the conditional distribution g(e\,p\\v\,0) reduces to 
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Figure 3. The 20 and 68 per cent percentiles of g(ei,pi \i>i, 0) for two different halo mass Mi = 2.2 X 10 12 (left panel) and Mi = 
3.4 X 10 13 M®/h (right panel). These mass scales correspond to a smoothing radius Ri = 2 and 5 h 1 Mpc, respectively. On large 
scale, the shear is smoothed on a fixed radius Rq = 10 /i _1 Mpc (Mo = 2.7 X 10 14 Mq/A). The shape of the large scale region is either 
a proto-cluster, the precursor of a filament or a sheet-like structure, characterised by (eo,poi^o) = (0.3,-0.1,1), (1.5,-1.3,0.5) and 
(1.3,1.1,0.4), respectively. The dashed lines indicate the boundary of the domain \pi\ > e\. The interior of the triangle bounded by 
(ei,pi) = (0, 0), (i, — j) and (i, 1) shows the region where A3 > 0. For a given mass Mi, protohaloes which collapse within filament-like 
or sheet-like structures are, on average, initially more asymmetric than those which will form in spherical-like environment. The average 
asphericity and the scatter increase with decreasing halo mass. 



the unconditional distribution, eq. @, derived by Doroshke- 
vich. It vanishes outside the domain defined by |pi| < e\. We 
stress that these expressions are valid for any window func- 
tion. In the particular case of a sharp fe-space filter however, 
the spectral parameter is simply 7 = 00 /cti . 

The distribution (|27[1 is sufficient to estimate the magni- 
tude of the environmental effect which arises from the statis- 
tics of the initial shear field. 

3.2 Statistical correlation between protohalo and 
environment in Gaussian initial conditions 

We illustrate how the local characteristics of the shear tensor 
smoothed on the scale of collapsing haloes depend on the 
initial geometry of their large scale environment. 

3.2.1 Shear ellipticity and prolateness 

First, it is worthwhile studying how the conditional prob- 
ability distribution changes with the correlation strength 
7. To this purpose, we evaluate the probability (|27[) f° r a 
density 5i = 2, a value sufficiently large so that the small 
scale overdensity collapses at moderate to high redshift re- 
gardless of the asymmetry parameters. For a filtering scale 
Rx = 2 /i _1 Mpc, i.e. a halo mass Mi ~ 2 X 10 12 M Q /h, 
this corresponds to a density contrast v\ = 1.28. On large 
scale Ro > Ri, we take the shear ellipticity, prolateness 
and density to be (eo,po, vq) = (1.4, —0.8, 0.6). These values 
are appropriate for a filament-like configuration. In Fig. [2] 
the 20 and 68 percentiles of the conditional probability dis- 
tribution g{e\,p\\v\,0) are plotted for 7 = 0, 0.4, 0.6 and 



0.8 (contours from top to bottom). A larger value of 7 in- 
creases the asymmetry of the distribution and sharpens its 
maximum. Recall that, at fixed value of R\, the correla- 
tion strength 7 decreases with increasing Ro (see Fig. 
For instance, with R\ = 2 h~ Mpc and a reasonable en- 
vironment size Ro ~ 10 h~ 1 Mpc, we have 7 ~ 0.6. The 
resultant distribution is displayed third from the top. In 
this case, the most probable values of the shear ellipticity 
and prolateness are (e m ,p m ) = (0.48,-0.15), significantly 
different than those in the limit 7^0, namely (0.35,0). 
Note that the top-hat smoothing artificially reduces the as- 
phericity of the large scale environment. Therefore, the val- 
ues of e m and p m certainly underestimate the asymmetry 
that could be measured from N-body simulations (Bardeen 
et al. 1986). Finally, unless otherwise stated, we shall adopt 
Ro = 10 fc _1 Mpc (Mo = 2.7 x 10 14 M e /h) as the environ- 
ment radius in the remaining of this paper. 

Fig. [3] further illustrates the environmental dependence 
which arises from the statistical properties of the shear ten- 
sor. Contours are plotted for three different configuration 
shapes of the large scale environment. The halo mass is now 
Mi = 2.2xl0 12 (left panel) and M a = 3.4xl0 13 M Q /h (right 
panel), and the resultant correlation coefficient is 7 = 0.58 
and 0.85, respectively. Also, we choose v\ so that 5i = 2 in 
all cases. The configuration shape of the large scale region 
is either a proto-cluster of signature (+ + +), the precursor 

of a filament (+ -j ), or a sheet-like structure (H ). 

Clearly, at fixed halo mass Mi, the protohaloes that col- 
lapse within the pancake or the filament are initially more 
asymmetric than those that will form in the spherical-like 
region. Furthermore, the asymmetry increases noticeably 
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with decreasing halo mass. For Mi = 2.2 x 10 12 M //i, 
the most probable values of ei and pi are (e m ,p m ) = 
(0.32, 0), (0.48, -0.14) and (0.39, 0.05) respectively, whereas, 
for Mi = 3.4x 10 13 M /7i, we find (e m ,p m ) = (0.17, -0.02), 
(0.37, -0.22) and (0.28, 0.14). In the high mass halo case, the 
conditional distributions exhibit much less scatter as a result 
of the larger values of 7 and v\. 

Since the density contrast is independent of the shear 
ellipticity and prolateness, the probability of finding a virial- 
ized halo in a given environment of initial shape (eo,po, va) 
is modulated by g{vi\vo), which is very sensitive to 7. For 
the low mass halo considered here, this probability is 0.38, 
0.31 and 0.17 for the cluster, filament and pancake config- 
urations, respectively. For the massive protohalo however, 
g(v\\vo) is significantly nonzero (0.10) for the spherical-like 
environment solely. Thus, low mass haloes form in the mildly 
overdense structures of the primeval density field, while the 
most massive collapse mainly in the densest, weakly aspheri- 
cal regions. This is a distinctive feature of hierarchical forma- 
tion models (e.g. Kaiser 1984; Mo & White 1996) that adds 
to decrease the scatter in the relation between the environ- 
ment and the local properties of the shear with increasing 
halo mass. This statistical effect may provide an explana- 
tion for the strong environmental dependence of low mass 
haloes. 



3.2.2 Alignment of principal axis frames 

It is fairly straightforward to derive probability distributions 
for the relative orientation of the principal axis frame from 
the results of i]3.1l Appendix [B] provides the details of the 
calculation. Briefly, we parametrise the rotation matrix R in 
terms of the Euler angles < ip,ip < 2n and < i9 < tt. We 
adopt the ZXZ convention so that $ is the angle between the 
two minor (third) axes. The trace tr (RAR T a) is then de- 
composed into a sum of rotation matrices Tj m m2 (ip,-d,ip). 
The integral over the variable ip can be performed easily and 
yields the conditional probability P(#, tp\0, 1) given the val- 
ues of eo, ei etc. In Fig. [3] contours of constant P($, <^|0, 1) 
are plotted for the cluster, filament and pancake-like config- 
urations described above. They are shown for the halo mass 
Mi = 2.2 x 10 12 (left panels) and Mi = 3.4 x 10 13 M©/ft 
(right panels) considered in Fig. [3] In all cases, the asymme- 
try parameters e\ and p\ assume their most probable val- 
ues. Unsurprisingly, the alignment between the minor axes 
of the shear smoothed on scale Ro and Ri is weaker for the 
low mass halo as a result of the smaller values of 7 and v\. 
At fixed halo mass however, the alignment is substantially 
stronger in the case of the filamentary structure (middle 
panels). This mostly follows from the fact that, for a con- 
figuration shape with one or two positive eigenvalues, the 
alignment is strongest along the axis of symmetry, which 
coincides with the minor axis for a filamentary structure. 
Conversely, the alignment between the major axes is much 
stronger in the pancake configuration. These results show 
that the principal axis frame of the tidal tensor smoothed 
on the protohalo mass scale cannot be assumed independent 
of that defined by the environment. This is especially true 
when the large scale configuration shape is highly asymmet- 
ric. Accounting for this correlation has a large impact on the 
probability distribution of the spin parameter (see £14.4. 3D . 




1 0.8 0.6 0.4 0.2 1 0.8 0.6 0.4 0.2 
cos(iS) cos(iS) 



Figure 4. Contours of constant probability P($, iy?|0, 1) for the 
cluster, filament and pancake-like configurations described in the 
text. They are shown for a halo mass M\ = 2.2 X 10 12 (left panels) 
and 3.4 X 10 13 Mq/Ii (right panels). The asymmetry parameters 
ei and pi assume their most probable value. Levels of contours 
decrease by a factor of 2. 1} is the angle between the minor axes. In 
all cases, the probability is highest along the vertical line cos 1? = 
1. The alignment is strongest in the case of the filament, whose 
minor axis coincides with the axis of symmetry. 



To summarise, a correlation between the local proper- 
ties of the shear and the configuration shape of the environ- 
ment is expected in Gaussian initial conditions because fluc- 
tuations on different scales are correlated. Our results nicely 
demonstrate the large magnitude of this statistical effect. At 
fixed halo mass, the shear of a perturbation that collapse 
within filaments or pancakes is on average more asymmet- 
ric than in spherical regions. The principal axis frame tends 
to be aligned along the axis of symmetry of the external 
mass distribution. Furthermore, the scatter in the relation 
between the primeval shear field and the geometry of the 
large scale environment increases strongly with decreasing 
halo mass. This will almost surely have a significant impact 
on the properties of virialized haloes in the ellipsoidal col- 
lapse model since the critical density threshold S cc depends 
strongly on the initial values of e\ and p\. We will further 
discuss the importance of this statistical correlation in the 
next Section. 
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4 ENVIRONMENTAL DEPENDENCE IN THE 
ELLIPSOIDAL COLLAPSE DYNAMICS 

We now concentrate on the dynamical origin of the envi- 
ronmental dependence in the ellipsoidal collapse model. We 
employ a simplified model based on the collapse of homo- 
geneous ellipsoids that takes into account the interaction 
between a collapsing halo and its surroundings. To antici- 
pate the results of this Section, we find that haloes residing 
in large overdensities virialize earlier. The environment den- 
sity is the key parameter in determining the virialization 
redshift. We incorporate this dynamical effect into the ex- 
cursion set formalism by means of a collapse barrier whose 
height depends on the environment density. This approach 
greatly simplifies the calculation of the environmental de- 
pendence of halo properties. It also predicts a clear correla- 
tion between formation redshift, large scale bias and envi- 
ronment. 



4.1 Homogeneous ellipsoidal dynamics 

4-1.1 Equation of motion 

We consider the collapse of a triaxial perturbation embedded 
in a uniform background assumed to be a ACDM cosmol- 
ogy. We neglect the influence of nonlinear substructures on 
the gravitational evolution. The initial (Lagrangian) volume 
occupied by this overdense fluctuation is a (uniform) sphere 
of comoving radius Ri. 

Following Peebles (1980) and Eisenstein & Loeb (1995), 
the (proper) position of any point interior to the ellipsoid is 
conveniently described as r a — A a/3 q 13 , where q T q < 1 and 
repeated indices are summed. The matrix A is a function 
of time solely. At all time, the equation defining the outer 
shell of the ellipsoid is r T (AA T ) r = 1. The principal 
axis lengths {Ak,k — 1,2,3} and directions of the ellipsoid 
can be found by diagonalizing AA T = QAQ T , where Q is 
orthogonal and A is a positive definite diagonal matrix. In 
this model, the potential is quadratic in the coordinates, 
$(r) = 1/2 $ a/3 r a r /3 , so that the external and internal forces 
preserve the homogeneity of the protohalo at all time. Intro- 
ducing the time variable r = ln(a) instead of t (e.g. Barrow 
& Silk 1981; Nusser & Colberg 1998), the equation of motion 
reads 

A Q/3 - (1 + ,(r))A^ = -~fi m (r) $ Q7 A 7/3 , (30) 

7 

and is manifestly independent of the Hubble constant. Up- 
per dots denote derivatives with respect to r, q(r) = 
fim(7")/2 — £2a(t) is the deceleration parameter (q(r) — 1/2 
in a EdS Universe) and the gravitational potential <E> is in 
unit of 4nGp m . Q m (z) and Sl\(z) are the matter and vac- 
uum density in unit of the critical density, respectively. The 
initial conditions are set by the Zeldovich approximation 
(see below). Virialization occurs when the three axes have 
collapsed. To prevent axis k from shrinking to arbitrary 
small sizes, we halt its collapse when Ak/aRi — f r = 0.177. 
This freeze-out radius is chosen so that the virialized object 
is 178 times denser than the background (Bond & Myers 
1996). When axis k has collapsed, we set the radial com- 
ponent of the velocity and acceleration in that direction to 
Ak — Ak = f r aR\. This way we end the radial collapse but 



leave the tangential velocity unchanged, so that the angular 
momentum of the protohalo is conserved (see, e.g., Eisen- 
stein & Loeb 1995 for details). 

4-1.2 External shear field 

The external force exerted by the rest of the Universe on the 
triaxial perturbation may be generically written in terms of 
the Green function G(r, r') = \r—r'\~ 1 .It is natural to adopt 
spherical coordinates as perturbations grow from an initially 
homogeneous and isotropic background. The potential inte- 
gral can thus be expanded as a multipole series (Binney & 
Tremaine 1987), 

*(r) = (21 + I)' 1 r l qf m {r) If* (ft) , (31) 

/ m 

where r= rh, Y ; m (n) are the spherical harmonic functions, 
and the coefficients 

qL(r) = /dV YF(h' )W) (r'y 1 - 1 (32) 

are the multipole moments that characterise the potential 
at any (interior) point r. The constant and the dipole 1 = 1 
term, which corresponds to a translation, do not alter the 
shape of the whole region and can be dropped out. The 
quadrupole / = 2 term describes the force distorting the cen- 
tral region. Analytic calculation (Quinn & Binney 1992) in- 
dicate that it dominates the higher order terms ignored here. 
This justifies to some extent the assumption of a quadratic 
potential. 

The quadrupole q£ m of the external shear will generally 
not be aligned with that of the protohalo region (see Sj3| . Al- 
though one could treat the nonlinear evolution of the exter- 
nal mass distribution with concentric shells (Chandrasekhar 
1969; Binney & Tremaine 1987; Ryden & Gunn 1987; Eisen- 
stein & Loeb 1995), we adopt a simpler approach and as- 
sume that the protohalo is embedded in a single, large scale 
triaxial region of initial comoving radius Rq > Ri, not nec- 
essarily lined up with the protohalo region. In principle, the 
protohalo will cause the large scale perturbation to warp, 
producing in return a non-quadratic potential that breaks 
the homogeneity of the small-scale ellipsoid. To avoid this 
problem, we assume that the large scale triaxial perturba- 
tion and the background are unaffected by the collapse of the 
protohalo and remain homogeneous at all time (Icke 1973; 
White & Silk 1979; Eisenstein & Loeb 1995). We can there- 
fore evolve the large scale ellipsoid independently using the 
model of Bond & Myers (1996). 

4-1.3 Gravitational potential of the protohalo 

In the principal axis frame of the large scale ellipsoid, we 
write down the total gravitational potential of the proto- 
halo as 

W — FRW + ^E,0 + *E,1 + ^zel \ M ) 

where $p£ w = (1 - 2n A (r)/O ro (r)) /3 8 a0 is the contribu- 
tion of the smooth background and 

•^o = b , a 
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Figure 5. Collapse redshift z c for the protohalo axes as a function 
of the initial ellipticity ei. We have assumed a ACDM Universe. 
The panels show how z c varies with the initial protohalo overden- 
sity 5i and the large scale density So and shear ellipticity eo. In 
all panels, po = Pi = and the dotted curve is the reference case 
(eo,5o) = (0,0). The solid and dashed curves show the collapse 
redshifts for (eo,5o) = (0, ±1) and (eo,<5o) = (0.3, ±1), respec- 
tively. The dotted-dashed curve also has (eo,5o) = (0.3, ±1), but 
the initial alignment of the shear is weaker (see text). 

* Z = 1^- <*°) £ b ^ Q Q7 Q' 37 (34) 

7 

are the gravitational potential associated with the large scale 
triaxial perturbation and with the remaining mass in the 
protohalo, respectively. The matrix elements Q Q! ' 3 (r) depend 
generally on r because the protohalo may be rotating. The 
functions b a (r) are defined in Appendix iJC] which provides 
details on the potential of an homogeneous, triaxial ellipsoid. 
It is worth emphasising that, since b a {r) are homogeneous of 
degree zero, the potentials (|34[1 do not depend on the value 
of Ro and Ri. So(t) and <5i(t) are the relative overdensity 
of the large scale environment and the protohalo, respec- 
tively. Since the contribution of the traceless part of both 
and is of second order only, we include the linear 
approximation for the external shear field, 

7 

where A zc i i7 = D(t) (A 7 — <5i/3) is a linear function of the 
initial eigenvalues A 7 of the shear tensor smoothed on scale 
R\. Qf = Q a,3 (Ti) describes the initial orientation of the 
protohalo relative to that of the large scale region. This 
ensures that the evolution consistently reduces to the Zel- 
dovich approximation at early times (Bond & Myers 1996). 
The initial conditions are explicitely 

A al3 (r t ) = ai Ri (1 - DiX ) Qf 

A Q/3 (rO = A Q V0 - a i R 1 D i \ l 3 Qf , (36) 

where at — a(r^) and Di — D(ri). Note that the initial 
tangential velocities are zero. 



The simplifications of this model notwithstanding, our 
calculation should provide a quantitatively useful descrip- 
tion of the impact of environment on the collapse of dense, 
triaxial regions. 



4.2 Effect of environment on the redshift and 
critical density for collapse 

Our first task is to study the dynamical effect of the large 
scale environment on the collapse of the protohalo. Once 
the cosmological background is chosen, the evolution of the 
protohalo is governed by the initial values of (ei,pi,Si), 
(eo,po,<5o) and the Euler angles (</?, i9, ip), which describe 
its initial orientation relative to that of the large scale el- 
lipsoid. For simplification, we choose po — pi — and set 
(<p,tf,ip) = (0.2,0.2,0.2) (in radian). Note that this partic- 
ular orientation has a reasonable probability of occurring 
regardless of the configuration shape of the environment 
(see Fig. 2} . The starting redshift is taken to be a hundred. 
We store the collapse redshift z c of the three well 
as the (linear) critical density S cc at virialization. Results 
are presented for a ACDM Universe with fl m — 0.238 and 
= 0.762. 

Fig. [5] examines how the collapse redshifts z c of the 
three protohalo axes change with the model parameters. 
For reasons that should become apparent below, z c is plot- 
ted against the shear ellipticity e\. In all panels, the dotted 
curve is the reference case (eo,5o) = (0,0). The solid and 
dashed curves indicate the collapse redshifts for (eo,<5o) = 
(0, ±1) and (eo,<5o) = (0.3, ±1), respectively. To highlight 
the effect of changing the relative orientation, we have also 
plotted as dotted-dashed curves the collapse redshifts for 
(eo,5o) = (0.3, ±1) and a weaker initial shear alignment, 
(ifij'&jip) — (1.5,1.5,1.5). Also shown are the values of 8q 
and 8i, linearly extrapolated to present epoch. For positive 
values of So , the tidal force exerted by the large scale pertur- 
bation delays the collapse along the first axis but enhances it 
along the third, reducing thereby the anisotropy that arises 
from the linear term, eq. (|35[) . By contrast, the difference 
between the collapse redshift of the major and minor axis 
is enhanced for So < 0. In all cases however, the collapse 
redshift of the intermediate axis is barely affected and re- 
mains close to the value predicted by the spherical collapse 
model (see Appendix of Shen et al. 2006). Clearly, at fixed 
overdensity Si, the haloes embedded in the large density 
environment virialize earlier. The strength of this effect in- 
creases with the shear ellipticity e\. It is very sensitive to 
the initial conditions. For So ~ — 1, an external quadrupole 
shear eo > delays the halo virialization as compared to the 
spherical case. This delay is still present for So = +1 when 
the halo collapse is initially close to spherical (The delay 
is most significant for the low value of Si), but the trend 
reverses when the ellipticity gets larger than e\ J>0.1. In ad- 
dition, changing the initial alignment can also increase or 
decrease the difference in collapse redshift. Notice also that, 
at low redshift, the increasing contribution of the cosmolog- 
ical constant to the energy density slows down the collapse 
of the third axis noticeably. 

When all the parameters but S\ are held fixed, then 
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Figure 6. Critical density 8 cc (ei,z) for the collapse of a proto- 
halo perturbation at z = in unit of S BC = 1.673, the critical 
density for a spherical collapse in the ACDM cosmology consid- 
ered here. <5 CC is plotted as a function of the shear ellipticity e\. 
The various curves indicate how 8 CC varies for the initial config- 
urations considered in Fig. [5] 'WA' labels the configurations for 
which the initial orientation of the shear principal axis frame is 
given by (ip, if)) = (1.5,1.5,1.5) (see text). The dotted curve 
shows the reference case (eg, So) = (0,0). Note that po = pi = 
in all cases. 



there is a unique value of Si = 5 ec (ei,z) which leads to 
the collapse of the minor axis (virialization) at redshift z. 
In Fig. [6] the critical density for collapse at z = is plotted 
in unit of S BC for the initial collapse configurations consid- 
ered in Fig. [5] We have also set po = pi = 0. As before, 
the dotted curve shows the reference case (eo,<5o) = (0,0). 
The label 'WA' designates the configurations that have a 
weaker initial alignment, (ip,-&,tp) = (1.5,1.5,1.5). Clearly, 
the critical density 5 ec is lower for larger values of <5o- This 
owes to the fact that, at fixed initial density Si, haloes that 
reside in relatively high density environments collapse ear- 
lier. Furthermore, Fig. [6] also indicates that the environment 
density So is the most influential parameter. The large scale 
ellipticity eo, for instance, has a significant impact on the 
critical density <5 CC only when So !>0. These results suggest 
that an average critical density B(ei, z, So) should provide a 
good description of this environmental effect if the scatter 
around it is not too large. We investigate this possibility in 
the remainder of this Section. 



4.3 Environmental dependence as moving barriers 

As seen in 331 the local properties of the shear depend sub- 
stantially on the large scale environment. Here, we consider 



1 In what follows, we shall omit writing the other variables, but 
recall that 5 ec is generally a function of the 12 variables that 
parametrise the shear eigenvalues and the relative orientation of 
the principal axis frames on scale Ro and R%. 



a large ensemble of initial halo-environment configurations 
and examine the resultant distribution of critical density 
S e c- We find that it is a reasonable approximation to use an 
average critical density B(ei, z, So) and neglect the scatter 
around it. We provide a fitting formula to B(ei, z, So) which 
facilitates the inclusion of the environmental dependence of 
the kind considered here into the excursion set formalism. 



4-3.1 Monte-Carlo simulations 

Chiueh & Lee (2001) have shown that random realisations 
of the linear deformation tensor can be simulated by draw- 
ing six independent Gaussian variables. Their algorithm can 
be easily extended to generate joint realisations of the shear 
tensor which satisfy the correlation property (JSJ . In the basis 
defined in eq. (|12[1 , the 12 x 12 covariance matrix M decom- 
poses into a direct sum of 2 x 2 block-diagonal matrices. As 
a result, the variables X = {u, v, w, £4, £5, ^} can be simu- 
lated using the following linear transformation 



(37) 



where, again, the subscripts and 1 indicate that the shear 
is smoothed (with a tophat filter) on scale Ro and Ri re- 
spectively. y + and y_ are two Gaussian random deviate of 
dispersion unity, a x = (X 2 ) is the rms variance in the 
variable X. For instance, a x — 1/VT5 when X = v. The 
remaining components of the shear tensor, £1, £2 and £3 are 
readily obtained from the linear relations (|12|l . 

Since the small-scale perturbation is identified as a halo 
at some redshift 2 > 0, while the exterior ellipsoid is assumed 
uncollapsed at that redshift, we constrain So so that it does 
not exceed 0.9 Si. We also enforce the constraint Si > 1.6. 
For a radius Ri — 2 /i _1 Mpc, this corresponds to a density 
threshold vi ;>1. In other words, the vast majority of our 
haloes form out of one standard deviation fluctuations. In 
practice, we generate random realizations of the initial con- 
ditions and reject the cases that do not satisfy these con- 
straints. 



4-3.2 Collapse barriers 

We generate a large ensemble of initial collapse configura- 
tions. We evolve each realisation separately using the model 
described in ij4.1l and store the value of the density contrast 
Si = S cc which corresponds to collapse at redshift z. 

The distribution of critical density S cc /S sc and shear 
ellipticity ei is plotted in Fig. [7] for several values of the 
primeval environment density, evenly spaced in the range 
— 2 < So < 1. Collapse occurs at z = 0. The crosses indicate 
the actual values of collapse densities and shear elliptici- 
ties of (a small subset of the) individual realisations. Also 
shown as the solid curve is our approximation to the criti- 
cal collapse boundary B = S cc (ei, z) defined by the implicit 
equation (Sheth, Mo & Tormen 2001) 



5 ec (ei, 



S BC (z) 



l + 0i 



5 ej 



5 e c(ei,z) 2 



S sc (z)i 



02 



(38) 
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Figure 7. Distribution of critical densities <5 CC and shear ellipticities e\ as a function of the environment density <5o for a collapse redshift 
z = 0. The crosses indicate the actual values of <5 ec and ei of 10 4 individual realisations. The solid curves show our approximation, 
eq. H38II . to the critical collapse boundary. The dashed curves indicate the (approximate) boundaries of the domain (ei,<5 ec ) sampled by 
the random realisations. They are the barrier shape l|38|l with 5q = ±5. 



The functions f3\ and /3 2 generally depend on both the red- 
shift 2 and the environment density So- We have adopted 
the simple functional form 

/3i(a,S ) = bi(l + z)- di exp(-aS ) • (39) 

The exponential factor ensures that f}\ and $2 are strictly 
positive. In addition, we have enforced the constraint di > 
to account for the (slight) decrease of the environmental de- 
pendence with redshift. Note also that, in the limit of large 
environment densities, the moving barrier B(e\, z) tends to- 
wards the constant spherical barrier B = S sc (z). The coef- 
ficients bi, a and di are found by fitting the barrier shape 
()38[) to the critical collapse densities of 3 x 10 5 realisations 
with environment density in the range — 2 < So < 1. We find 

61 « 0.412, ci w 0.113, di » 0.0576 

b 2 « 0.618, c 2 « 0.0451, d 2 « 0.0485 . (40) 



For So = and a = 1, we obtain /3i = 0.412 and /3 2 = 0.618, 
in good agreement with the values inferred by Sheth, Mo & 
Tormen (J3i = 0.47, /3 2 = 0.615). A visual inspection of Fig. □ 
has convinced us that the barrier shape (|38[l provides a good 
approximation to the increase of the critical collapse density 
with ellipticity and to its dependence on the environment 
density in the range < ei < 0.45. To guide the eye, we 
have also plotted as dashed curves the (approximate) upper 
and lower boundaries of the domain (ei, <5 CC ) sampled by the 
random realisations. These boundaries are the barrier shape 
(1381) with an environment density 80 = ±5. 

It should be noted that the barrier (|38p is in good 
agreement with exact result for po = Pi — 0. In this case, 
the parametrization of Sandvik et al. (2007) and Sheth et 
al. (2001) agree fairly well, although the former describes 
the two-dimensional surface 5 ec (ei,pi) more accurately. For 
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Figure 8. Top panel : Distribution of critical density and shear 
ellipticity for three different protohalo radii i?i = 1 (circle), 2 
(square) and 5 h~ 1 Mpc (triangle). The collapse redshift is z = 
and the environment density is <5o = 0. The average density for 
collapse, eq. (138 b . is plotted as the solid curve. Bottom panel : 
Distribution of initial shear ellipticity and prolateness for the in- 
dividual realizations shown in the top panel. The interior of the 
triangle shows the region where A3 > 0. The big empty symbols 
indicate the peak of the distributions. 



So 7^ however, we found that the critical density for col- 
lapse depends on the shape parameters in a complex way. 
This is the reason why we resort to the Sheth et al. func- 
tional form. 

Following Sheth et al. (2001), we shall interpret (f38)l as 
a 'moving' barrier. Such an interpretation has the advan- 
tage that, once the barrier shape is known, the excursion 
set formalism can be use to quantify the dependence of halo 
properties on environment. Moreover, it is computationally 
more efficient that studying the first crossing distribution of 
multi-dimensional random walks (e.g. Chiueh & Lee 2001; 
Sheth & Tormen 2002; Sandvik et al. 2007). 



4-3.3 Mass scale-ellipticity relation 

Before we examine the impact of this dynamical interaction 
on the properties of collapsed haloes, we need to express the 
critical density (|38fl in terms of the halo mass M\. Thus far, 
we have only considered the collapse of regions with initial 



radius R\ — 2 /i _1 Mpc. The top panel of Fig. [8] shows the 
distribution of critical densities o" ec for three different val- 
ues of Ri. The average density for collapse is plotted as the 
solid curve. Results are shown for an environment density 
So = only. Note, however, that we have repeated this cal- 
culation for other values of So and found good agreement 
between the critical densities of individual realisations and 
the approximation (|38[) . This confirms that most of the envi- 
ronmental effect seen here arises from variation in the large 
scale density So- At fixed Ro, changing R± merely affects the 
scatter around the collapse boundary B(ei,z), unsurprising 
since the functions 6t(r) that characterise the potential of 
the protohalo and its environment are independent of Ro 
and R\. The decrease in scatter with increasing Ri is a di- 
rect consequence of the statistical correlation explored in 
[J3] This is clearly seen in the bottom panel of Fig. [8] where 
the distribution of initial shear ellipticity and prolateness 
is plotted as a function of Ri. The filled symbols indicate 
the actual, most probable values (e m ,p m ). They increase 
(monotonically) with the protohalo radius Ri . This suggests 
relating the mass scale Mi to the expectation value of the 
asymmetry parameters. 

In line with the interpretation of Sheth et al. (2001; see 
also Shen et al. 2006; Sandvik et al. 2007), we use the av- 
erage values ([TJl to translate the product eiSi into a peak 
height v(Mi,z) = <5 sc («)/cr(Mi) B Gao & White (2007) 
have shown that the properties of haloes in the Millennium 
simulation obey this scaling relation over a large redshift 
range. In terms of this scaled variable, the collapse bound- 
ary eq. (|38p becomes 



B(u,z)=S ec (z) 1 + 0! 



-/3 2 



(41) 



Sheth & Tormen (2002) provide an analytic fit for the first 
crossing distribution associated to this barrier shape, which 
allows us to easily calculate halo properties such as forma- 
tion redshift and bias. 



4.4 Environmental dependence of halo properties 

The main focus is to quantify the environmental dependence 
of halo properties arising from the moving barrier, eq. (|38p . 
Even though the ellipsoidal collapse model does not provide 
an excellent description of the simulations, it does neverthe- 
less a much better job than the spherical collapse (Sheth & 
Tormen 1999; 2002). Therefore, the validity of a compari- 
son between predicted and observed environmental effects 
should be preserved. We calculate the distribution halo of 
halo formation redshift and large scale bias associated to 
that collapse boundary. We end this Section with a discus- 
sion of the halo spin parameter. 



4-4-1 Environment density and formation redshift 

The halo formation redshift z[ orm is commonly defined as the 
epoch at which the main progenitor has accumulated half of 



2 The peak height u is the typical amplitude of fluctuations that 
produce haloes of mass M\ by redshift z. A characteristic mass 
for clustering M*(z) can then be defined through v(Mi,z) = 1. 
For the present cosmology, M*(0) fa 2.6 X 10 12 Mq / h. 
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Figure 9. The differential probability distribution of formation 
redshift, P(zf orm ), as a function of halo mass and (Lagrangian) 
environment density So- P(2[ orm ) has been computed from the 
moving barrier l|41j l. The halo mass is given in unit of the char- 
acteristic mass Af* w 2.6 X 10 12 M Q /h. For M < M t , we also 
indicate Azf orm , the difference between the lowest and highest 
median formation redshift. The mean formation redshift 2f orm in- 
creases with decreasing So- The effect is strongest for small mass 
haloes M < M*. 



its final mass. According to the EPS theory, the probability 
that the formation redshift of present-day haloes of mass M 
(we will henceforth drop the subscript 1) is larger than z is 
given by (Lacey & Cole 1993) 



P(> Zf) 



s 2 



dS' 



M 



M'(S' 



■f(S'\S) 



(42) 



where £2 = S(M/2). When the excursion set theory is com- 
bined to the ellipsoidal collapse dynamics, the conditional 
first crossing distribution f(S'\S) shall be replaced by the 
following analytic formula, 



f(S'\S) 



where 



T(S'\S) = 



\T(S'\S)\ 



2n(S' - S) 3 /- 



■ exp 



(B(S') - B(S)f 



2(8' - S) 



(S-S') n d n [B(S')-B(S)] 
n\ dS' n 



(43) 



(44) 



This Taylor expansion provides a good fit to the conditional 
up-crossing probability for moving barriers of the form (|41[) 
(Sheth & Tormen 2002). Note that the ellipsoidal collapse 
model predicts too many haloes with high formation red- 
shift. The distributions are also broader than seen in the 
simulations (Lin, Jing & Lin 2003; Giocoli et al. 2007). 

The differential probability distribution P(Z{ otm ) = 
dP(> 2form) /deform is plotted in Fig. [9] for various halo 
mass M and (Lagrangian) environment density So- Clearly, 
the dependence of halo formation redshift on environment 
increases with decreasing halo mass. Equation (|41fl indeed 
shows that, for haloes of mass M ;>M*, the critical den- 



sity for collapse is B(v,z) w S sc (z), weakly dependent on 
environment density. This is the reason why the forma- 
tion redshift of massive haloes is nearly insensitive to So- 
By contrast, random walks associated to small mass haloes 
M <C M* cross the collapse boundary B(v,z) at relatively 
larger values of v and thus induce a stronger environment 
effect. For M = 0.01M+, we find a median formation red- 
shift of Zform = 1.19 for an environment density So = 1. This 
should be compared to Zf orm = 1-42 when So = —2. Notice 
that the probability distribution P(zf orm ) is more sensitive 
to the exponent P2 than the multiplicative factor f3\. The 
former contributes about two thirds of the environmental 
effect (see £|5.1l for more details). 

To allow for a direct comparison of our results with the 
analyses of N-body simulations, we need to express the envi- 
ronmental dependence of halo formation redshift in terms of 
the evolved Eulerian density. Mo & White (1996) and Sheth 
(1998) have shown how this may be accomplished within the 
spherical collapse model. The spherical collapse dynamics 
provides a relation between the Lagrangian radius Ro and 
density 80 and their Eulerian counterparts R and 8. In this 
model, the mass interior to each perturbation is constant, 
giving Ro = R(l + 5) 1//3 if one assumes that the primeval 
density fluctuations are small. Furthermore, the linear den- 
sity 80 is a monotonically increasing function of the present 
overdensity fluctuation 5 only. Mo & White (1996) have ob- 
tained an accurate approximation to this relation for an EdS 
Universe, 



•So (A) = 



1.686 



1.686 



1.35 



1.124 0.788 



+ ■ 



^2/3 AV2 A' 



(45) 



where A = 1 + 5. This interpolation formula is also valid 
in the ACDM cosmology considered here provided that 5 is 
not too large (S <J10). Ideally, we should calculate the rela- 
tive numbers of patches with (Ro, So) that have now evolved 
into regions (R,S) (see Sheth 1998). We should also take 
into account the triaxiality of the large scale environment in 
the conversion of S into So- Henceforth however, we will ne- 
glect these complications and use the spherical approxima- 
tion to relate the Lagrangian density to the Eulerian density 
at z = 0. This is sensible since, as we have seen, the linear 
overdensity So is the key parameter governing the correlation 
between collapse densities and environment. For illustration, 
the Lagrangian density is in the range —3 <C<5o <4 when the 
Eulerian density varies between 0.3 < (1 + S) < 5. 

Upon these assumptions, we find that the median for- 
mation redshift changes by 



A0 form w 0.07 when M = M* 
Aziorm » 0.33 when M = 0.01M* 



(46) 



when the evolved density varies in the range —0.7 < 8 < 4. 
Again, the statistical correlation induces an effect greater for 
low mass haloes M <,M* while, for M 3> M*, the offset in 
median formation redshift is barely discernible. In addition, 
the difference is of the magnitude seen in N-body simula- 
tions, where Az{ OT m ~ 0.3 — 0.5 at M ~ 0.01M* for a large 
scale overdensity varying in that same range (Harker et al. 
2006). On the other hand, our model predicts that haloes 
in denser regions have a lower formation redshift than those 
in less dense regions. This is opposite to the behaviour re- 
ported by Harker et al. (2006), who find that, in overdense 
regions, haloes in the densest environment assemble earliest. 
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Figure 10. The large scale bias factor b(u) as a function of the 
peak height v = S BC /a(M). The shaded area indicates the am- 
plitude of the large scale bias when the Eulerian environment 
density varies in the range —0.7 < 8 < 4. The dashed curve is 
b(u) when <5 = 0, whereas the dotted line indicates the predic- 
tion of the spherical collapse. The upper and lower solid curves 
show the bias factor of haloes which have relatively high and low 
formation redshifts (see text for details). Haloes that assemble 
relatively early are more clustered and populate the less dense 
regions. 

A better treatment of the relation between Lagrangian and 
Eulerian regions should not affect these conclusions. 

Interestingly, however, these authors find that, in the 
most underdense regions S <^ — 0.4, the average formation 
redshift increases with decreasing environment density. De- 
spite the small number of haloes and the large scatter in 
formation redshift, the trend is robust, yet smaller than the 
effects present in the high density regions (Geraint Harker, 
private communication). This may be a manifestation of the 
environmental dependence discussed here. 

4-4-2 The age-dependence of halo bias 

Having established a correlation between formation red- 
shift and environment density, we now turn to the age- 
dependence of the halo bias. 

Mo & White (1996) and Sheth & Tormen (1999c) have 
shown that there is a direct relation between the halo bias 
and the shape of the collapse barrier. For a barrier of the 
form l|4ip . the large scale bias in Eulerian space can be ap- 
proximated by (Sheth et al. 2001) 

6(1/) = 1 + ^- \v 2 +l3 1 v 2 - 2 ^ 

v 2h 1 
~u*to (1-/32) (1-^/2) J ■ (47) 

This bias relation is plotted in Fig. [10] as a function of the 
peak height v. The dashed curve shows the halo bias at 
mean Eulerian density 5 — 0, namely in the case of the 



ellipsoidal collapse of Sheth, Mo & Tormen (2001), while 
the dotted line is the prediction of the spherical collapse. 
The shaded area indicates the amplitude of b(u) when the 
environment density varies in the range —0.7 < S < 4. The 
bias is roughly 10 per cent larger for the haloes residing in 
the most underdense region S = —0.7. 

To compare our results as directly as possible with those 
of Gao, Springel & White (2005), it would be useful to es- 
timate the bias factors of haloes which lie in the upper and 
lower tail of the formation redshift distribution. In prac- 
tice, however, such a calculation proves difficult in the EPS 
formalism. Instead we will consider a simpler approach mo- 
tivated by the relation between barrier height and median 
formation redshift. We have so far neglected the presence 
of stochasticity in critical collapse densities. In this respect, 
Fig. [7] shows that, at fixed environment density, there is a 
fairly large scatter in 5 ec . To estimate the importance of 
this scatter in the bias, we will also compute b(u) for the 
barriers that bound the locus traced out by the random re- 
alisations when the environment density varies in the range 
—0.7 < S < 4 (i.e. a Lagrangian density —2 < So < 1). 
These upper and lower barriers are shown in Fig. [7] as the 
dashed curves. They induce a formation redshift distribu- 
tion strongly biased towards the extremes. For M = 0.1M*, 
the median formation redshift associated to the upper bar- 
rier is larger by Azf orm = 0.51 than the value obtained from 
the lower barrier. We use these barriers to define our 'old' 
(high Ztorm) and 'young' (low 2form) haloes. Of course, this is 
a crude approximation to the 10 (20) per cent tails consid- 
ered by Gao et al. (2005), but we have not found any better 
alternative. 

In Fig. 1101 we show the resultant bias relations as the 
solid curves labelled by 'low Zform' an d 'high Zf orm '. The rel- 
ative bias of our old versus young haloes increases smoothly 
with decreasing halo mass. The effect becomes large for v <A 
because of the considerably stronger dependence of the bar- 
rier shape on the peak height. In this regime, the large scale 
bias of the old haloes is ~50 per cent larger than for the 
young ones, an increase roughly twice as small as seen in 
the simulations. Overall, the behaviour is similar to that re- 
ported in Gao et al. (2005), where the correlation between 
halo clustering and formation history is strong for haloes less 
massive than M* only. A better modelling of the properties 
of haloes lying in the tails of the formation redshift distri- 
bution is required to make a more quantitative comparison 
between our predictions and the measurements of Gao et 
al. (2005). The point of the present analysis is to show that 
the moving barrier (|41[) could lead to a correlation similar 
to that seen in N-body simulations if the scatter in collapse 
density is taken into account. 

In contrast to the behaviour seen in N-body simulations, 
the predicted bias b(v) decreases with increasing environ- 
ment density, reflecting the flattening of the barrier shape 
for large values of So- In the limit j > 1, b(v) tends to- 
wards the value predicted by the spherical collapse. Since, 
on average, dense regions at the present time formed from 
relatively dense regions in the primeval fluctuation field, an 
anti-correlation between halo bias and large scale Eulerian 
density is expected in our model. This seems at first surpris- 
ing but, as recognised by Abbas & Sheth (2007), the large 
scale bias is not necessarily a monotonically increasing func- 
tion of environment density. Their results strongly suggest 
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Figure 11. Top panel : Effect of including the correlation be- 
tween the initial orientation of the protohalo and large scale en- 
vironment on the probability distribution of the spin parameter 
A. The left histogram shows P(X) when the correlations in the 
initial alignment are included, the right histogram when they are 
not. Bottom panels : Probability distribution of the cosine of the 
angle between the angular momentum of the collapsed halo and 
both filament direction (bottom left) and sheet normal vector 
(bottom right). A random distribution would be a flat line at 
P(cost?) = 1. Results are shown for haloes that collapse at z = 0. 
The histograms were drawn from ~ 10 4 random realizations of 
the initial conditions. 



that haloes in extremely underdense environment are more 
clustered than the mean. 



4-4-3 Alignment of halo spin parameter 

The dimensionless spin parameter A measures the amount 
of rigid rotation acquired by the triaxial perturbation before 
virialization. It is defined as (Peebles 1969) 



A 



W\E\ 
GM 5 / 2 



(48) 



where L = LI is the halo angular momentum, I is a unit vec- 
tor, and E — E pot + -Ekin, M are the total energy and mass 
of the ellipsoid, respectively. The quantities L and E can 
be expressed in terms of the matrix elements A a/3 (Peebles 
1980; Binney & Tremaine 1987; Eisenstein & Loeb 1995). 
In general, the energy E of the collapsing region is not con- 
served, notably because the kinetic energy of the ellipsoid 
is altered when an axis collapses. For simplicity, we use the 
last value of the energy before the first axis collapses in the 
calculation of A = jA| (Eisenstein & Loeb 2005). This is a 
good approximation as the change in total energy is usually 
small during the collapse. Note also that the resulting er- 
ror should be relatively small because A depends on yp5| 
solely. 

The top panel of Fig. [11] illustrates the effect of in- 
cluding the correlation between the initial orientation of the 
protohalo and its large scale environment on the probability 



distribution P(X) of the spin parameter. The left histogram 
shows P(X) when the correlation in the initial alignment is 
included, the right histogram when statistical independence 
is assumed. The distributions were drawn from 10 4 Monte- 
Carlo realizations of the initial conditions. Clearly, the in- 
clusion of correlations in the alignment of the principal axes 
has a strong impact on P(X) : It lowers the median spin by 
<;40 per cent. As pointed out by several authors (Steinmetz 
& Bartelmann 1994; Catelan & Theuns 1996; see also Hoff- 
man 1988), the assumption of random relative orientations 
overestimates the growth of angular momentum because, if 
the principal axis frames tend to be partially aligned, the 
angular momentum gain is reduced. Fig. [11] confirms that 
the correlation in the primeval alignment is an important 
factor, in agreement with the analysis of Lee & Pen (2001). 

Our median spin value is A mc( j ~ 0.005, an order of 
magnitude lower than those found in numerical simulations, 
where A mcd ~ 0.03-0.05 for haloes of mass M ^lO 11 M Q /h 
(Barnes & Efstathiou 1987; Bullock et al. 2001; Bett et al. 
2007). There are several reasons for this discrepancy. First, 
the number density of haloes of a given mass depends notice- 
ably on the environment density S : haloes are preferentially 
found in mildly overdense regions. Here, however, the distri- 
bution of evolved density 8 associated to the random reali- 
sations of Fig. [TT]is not representative of a fair halo sample : 
it peaks around S — 0, where the average spin at collapse 
is much lower than in high density regions S J>1. Second, 
the angular momentum of a collapsing region satisfies the 
equation 

kl-rlj 



Li 



-H- L e ijk T Kl V' 



(49) 



where T ij = (3/2) H 2 D. m (^ + &J cl ) is the torque and 
P J = (l/5)M(AA T ) ij is the inertia tensor. We note that 
the angular momentum vanishes to first order if the iner- 
tia tensor is zero at initial time, or if I and T are per- 
fectly aligned (White 1984; Catelan & Theuns 1996). In this 
limit, the growth of angular momentum is of second order 
only, L oc a 5 ^ 2 , and is dominated by nonlinear effects after 
turnaround (Peebles 1969). This is indeed the case here as 
we are considering perturbations that are initially spherical. 
A more realistic treatment should include first-order tidal 
torquing. Third, our axis collapse condition strongly reduces 
the amount of angular momentum gained by the halo after 
the second axis collapses. A different prescription, such as 
the one adopted by Eisenstein & Loeb (1995), would in- 
crease the magnitude of the spin parameter at third axis 
collapse. Finally, analytic calculations indicate that a signif- 
icant fraction of the angular momentum is acquired through 
the collapse of the outermost shells (Ryden 1988; Quinn & 
Binney 1992) . This can only be taken into account by a de- 
tailed modelling of the matter density and velocity profiles 
around the collapsing haloes. 

In spite of these limitations, it is worthwhile looking at 
the alignment between halo spin and environment as this 
quantity does not depend on the magnitude of the angu- 
lar momentum gained during the collapse. In the bottom 
panels of Fig. 1111 the histogram shows the alignment distri- 
bution at collapse time. The angle is measured between the 
angular momentum of the collapsing region and the sym- 
metry axis/plane of the mass distribution. Clearly, haloes 
show a strong tendency to have their spin aligned perpen- 
dicularly to the filament or parallel to the mass sheet. This 
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is in good agreement with the findings of Hahn et al. (2007) 
and Aragon-Calvo et al. (2007a), although the former au- 
thors found a clear correlation for haloes residing in sheet- 
like structures only and the latter a mass- dependent spin 
orientation in filaments (see also Sousbie et al. 2007; Tru- 
jillo, Carretero & Patiri 2006). More precisely, the alignment 
is strongest along the second principal axis of the shear ten- 
sor. This owes to the fact that, once the first axis has col- 
lapsed, Li ~ tijk o ^ ne S row th ra t e i s largest for the 
intermediate axis, L2 oc (qi — 0:3) I 23 , since the difference 
0.1 — «3 dominates the other two (see also Lee & Pen 2000). 
It is unclear whether the initial alignment between inertia 
and deformation tensors reported by Lee & Pen (2000) and 
Porciani, Dekel & Hoffman (2002b) can produce a similar 
correlation. We have not investigated this issue any further. 

To summarise, we have demonstrated that our simpli- 
fied model, which takes into account both the dynamical 
and statistical aspects of the ellipsoidal collapse, produces a 
clear correlation between formation redshift, large scale bias 
and environment density. The strength of the effect is of the 
same magnitude as seen in simulations. It is largest for low 
mass haloes, M <C M*, and fades as we go to high masses, 
M > Mi,. Haloes that formed at high redshift are substan- 
tially more clustered than those that assembled recently. 
This is precisely the behaviour reported by Gao, Springel 
& White (2005). On the other hand, this model predicts 
a negative correlation between formation redshift and envi- 
ronment density, in contradiction with the trend measured 
by Harker et al. (2006). However, simulations indicate that 
halo properties depend on environment in a complex way : 
in relatively underdense regions, the average formation red- 
shift increases with decreasing environment density (Harker 
et al. 2006) while haloes are more strongly clustered than 
the mean (Abbas & Sheth 2007). Our model produces an 
effect in the right sense. This suggests that the ellipsoidal 
collapse may apply in underdense regions where nonlinear 
effects are weak or absent. 



5 DISCUSSION 

In this Section, we discuss the sensitivity of the environ- 
mental effects considered above to the shape of the collapse 
barrier, as well as non-Markovianity and tidal interactions 
as potential sources of environmental dependence. 

5.1 Sensitivity to the barrier shape 

The strength of the environmental effect explored in H4.4I 
may somewhat depend on the parametrisation adopted for 
the moving barrier. While this is probably true for haloes 
of mass M J>M*, the detailed shape of the barrier should 
have a little impact when M <CM*. In this mass range, the 
average collapse density B(v, z), a Taylor series without loss 
of generality, should be dominated by the term of largest de- 
gree. For the simple case considered here, B(y, z) oc P\v~ 2 ^ 2 
when v 3> 1. 

3 Here the ct^s are the (ordered) eigenvalues of the large scale 
potential <E>g Q 
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Figure 12. Sensitivity of the median formation redshift to the 
shape parameters /3i and $2 ■ Contours of constant Zf orm are plot- 
ted for a halo mass M = 0.1M*. The contour levels are evenly 
spaced, with an interval Azf orm =0.1. The filled symbols indicate 
the median formation redshifts for the two extreme cases So = 1 
and —3 when fix and 02 assume the functional form l|39j l. The 
effect of varying /3i vanishes along a critical line P2 = /3S 1 where 
/3£ ~ 0.58 (dotted curve). 



The exponent $2 has a great influence on the correla- 
tion between formation redshift and environment density. 
This is clearly seen in Fig. 1121 where contours of constant 
•Zform are plotted for a halo mass M = 0.1M*. Shown for 
illustration are the median formation redshifts of two ex- 
treme cases So = 1 and —3 (filled symbol) when /3i and 
P2 assume the functional form (|39|l . The parameter f3i has a 
substantial impact on the formation redshift only when (52 is 
far from a critical (empirically determined) value /3§ — 0.58. 
Along this curve, the effect of varying /3i vanishes, presum- 
ably because the effective spectral index n e ff that controls 
the shape of M(S) in eq. (|42|l conspires to maintain the me- 
dian formation redshift constant. Since n c g varies with the 
mass scale, we expect (3% to change somewhat with the halo 
mass. The fact that the ellipsoidal dynamics predicts a value 
/?2 = 0.618 close to P2 is a coincidence. 

Regarding the environmental dependence of halo bias, 
note that, in the limits v 3> 1 and 1/ « 1, the bias offset 
is Ab(u) w -2/3 1 A/32lnuu 2 - 2 ' 32 /S sc to first order in Afc 
and A/32- In this regime, changing the value of P2 has a 
large effect on the bias of haloes lying at the extreme ends 
of the mass range. However, unlike the environmental de- 
pendence of formation redshift, it is the parameter (5\ that 
influences most the bias when the halo mass is in the range 
0.1 <,M/M+ <;10. We also note that the fitting formula (O 
of Sheth, Mo & Tormen (2001) holds for j3 2 strictly less than 
one only. It would be prudent to check again their approxi- 
mation, namely, derive the first crossing distribution from a 
large ensemble of random trajectories when the values of /3i , 
/?2 differ significantly from those at mean density (<5o = 0) . 

We have shown that each of the shape parameters /3i 
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and $2 has a distinct impact on the halo formation red- 
shift and large scale bias. Therefore, an environmental de- 
pendence of formation redshift and bias will be present re- 
gardless the exact values of these variables. 

5.2 Environmental effect from non-Markovianity 

In the excursion set formalism, a spherical symmetric win- 
dow function W(R, k) is used to define the trajectories of 
the linear density field 5(R) as a function of smoothing scale. 
When the sharp fc-space filter is adopted, S(R) executes a 
random walk. The property of Markovianity has been ex- 
ploited extensively to obtain analytic expressions for the first 
crossing distribution etc. This is the reason why we have con- 
sidered that particular window function in the computation 
of halo formation redshift and bias. However, the Markov 
nature of Brownian walks prevents any correlation between 
large scale environment and assembly history (see White 
1996 for a discussion). Halo merger trees are indeed non- 
Markovian across short time steps (Neistein & Dekel 2007). 
A natural way to introduce correlations would be to use an- 
other window function. This possibility has been considered 
by several authors (e.g. Peacock & Heavens 1990; Bond el 
al. 1991; White 1996; Schiicker et al. 2001; Nagashima 2001; 
Amosov & Schiicker 2004; Zentner 2007). Here we argue that 
non-Markovianess should lead to an effect that is larger for 
high mass haloes. 

For concreteness, let us consider the collapse of a per- 
turbation of comoving size Ri . To investigate the impact of 
environment on its formation history, we restrict ourselves 
to trajectories that obey the following two constraints : 

Co : S(R ) = uo 0o 

Ci : 8{R 1 ) = v 1 a 1 , (50) 

where, again, at — a(Ri) and Ro = 10 /i _1 Mpc is the scale 
of the large scale environment. We neglect triaxiality and 
choose v\ = 8 BC /<j\ to ensure that the halo has just collapsed 
by redshift z — 0. vo can be positive or negative, depending 
on whether the halo forms in a high or low density region. 

We can calculate the most probable trajectory S(R) 
given the constraints {Ci}. Since these are linear functional 
of the density field, the probability of possible realization 
S(R) can be expressed as a shifted Gaussian around an en- 
semble mean field (see Adler 1981; Hoffman & Ribak 1991; 
Van de Weygaert & Bertschinger 1996 for a rigorous treat- 
ment) 

8(R) = bW&c, , (51) 

where (i(R) = {6(R)d) is the cross-correlation between 
the field and the ith constraint, and £ij = (dCj) is the 
constraints' correlation matrix. The residual field 8 — 8 — 8 
is a Gaussian random field which is not homogeneous nor 
isotropic, but whose statistical properties are independent of 
the {Ci} (Hoffman & Ribak 1991). We define a normalised 
cross-correlation between the constraints and the field, 

1 f°° 

Ci(R)= / d\nkA 2 s {k)W(Ri,k)W(R,k) , (52) 

a ^ Jo 

where a = a(R) and W(R, k) is the tophat filter. The cal- 
culation of the matrix elements Qj is immediate. The mean 
field 5 can be expressed as 




R / tr'Mpc 



Figure 13. Mean field S(R) as a function of the smoothing ra- 
dius R for two different halo size Ri = 1 and 5 /i. —1 Mpc. The 
mean field obeys the following two constraints : i) 6(Ri) = S sc 
ii) <5(-Ro) = ±o"o. The solid and dashed curves show the aver- 
age trajectories which satisfy 8(Ro) = -co and 8(Ro) = +o"o, 
respectively. The shaded areas indicate the 68% scatter around 
the mean. For clarity, they are attached to the curves with 
5(Ro) = — oo only. The horizontal line is the constant spheri- 
cal collapse barrier 8 = 8 BC . Results are shown for the tophat 
filter. 

^<i^>MS-H('-^)}' <*> 

The variance of the residual is d~ 2 (R) = a 2 — Cij Ci an d 
does not depend on the constraints Co and C\. For a given 
mass scale Mi oc R\, the difference <5+ — <5_ between the 
mean field of two different large scale environments with 
density contrast v+ and V- < V+ is 

6+{R) - S-(R) = ~ O ° ^iffi ■ (54) 

This difference is negative (positive) when the smoothing 
radius is R < R\ (R > Ri). In other words, the density 
profile around 5(R 1 ) is steeper in low density regions. Note 
that, for a sharp fe-space filter, 5+ — <5_ vanishes when R < 
Ri since £i = Co/7, but is generally non-zero in the range 
R > Ri. 

Fig. [T3] shows the mean field S(R) for the tophat filter. 
8(R) is plotted as a function of the smoothing radius for two 
different halo mass Mi = 2.7 x 10 11 and 3.4 x 10 13 M e /h. 
On larger scale Ro = 10 h~ 1 Mpc, the density contrast is 
set to v+ = 1 (dashed curves) and V- = — 1 (solid curves). 
To ensure that the trajectories 8(R) describe the formation 
history of M = Mi haloes, we should also have constrained 
the density so that S(R) does not cross the barrier 8 = 8 SC 
on scale larger than Ri . Here we ignore this constraint since 
its implementation is not straightforward when the window 
function differs from a sharp fc-space filter. This is the rea- 
son why a substantial fraction of the trajectories penetrate 
the barrier 8 — S sc on scale R > R\. This caveat notwith- 
standing, the present calculation is adequate to understand, 
at least qualitatively, the effect of non-Markovianity. No- 
tice that the difference 8+ ~ 8- is approximately constant 
throughout most of the halo formation history (R/Ri < 1). 
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Most importantly, the effect is about twice as large for the 
high mass halo while, at fixed R/Ri, the scatter of the resid- 
ual field is lower by ~ 60 per cent. 

This suggests that the environmental dependence which 
arises from non-Markovianess is stronger for high mass 
haloes. This has also been pointed out by Zentner (2007). 
Therefore, it is unlikely to explain the trend seen in over- 
dense regions where, undoubtedly, correlations are stronger 
for low mass haloes (e.g. Gao & White 2007). However, it 
may apply for isolated haloes in voids or underdense regions. 
It should also leave a signature distinct from the ellipsoidal 
collapse which, as we have seen, induces a stronger depen- 
dence for low mass haloes. 



5.3 On tidal interactions 

Irregularities in the mass distribution induce non-radial mo- 
tions that slow down the collapse (Peebles & Growth 1976; 
Davis & Peebles 1977; Peebles 1990) . This "previrialisation" 
conjecture is supported by the numerical investigations of, 
e.g., Barrow & Silk (1981); Szalay & Silk (1983); Villumsen 
& Davis (1986); Lokas et al. (1996); and by the analytic cal- 
culations of Del Popolo et al. (1998; 2001), which indicate 
that tidal heating can counterbalance the effect of the shear 
and delay the collapse. Recently, Avila-Reese et al. (2005), 
Maulbetsch et al. (2006), Wang et al. (2006) and Diemand 
et al. (2007) have proposed that the assembly bias seen in 
N-body simulations originates from tidal interactions with 
a larger neighbour. They have shown that, at late time, the 
tidal field of massive neighbouring clumps halts the growth 
of haloes and, in many cases, even reduces their mass. Fur- 
thermore, tidal effects appears to have a larger influence on 
small mass haloes. Therefore, this could also explain why 
the age-dependence of clustering is stronger for low mass 
haloes. The large effect measured by Diemand et al. (2007) 
indicates that, in high density regions, tidal interactions are 
likely to overwhelm the environmental dependence arising 
from anisotropic collapse, and increase the average forma- 
tion redshift in overdense regions. 

Although the impact of tidal stripping can only be rigor- 
ously quantified with numerical simulations, the suppression 
of mass accretion through tidal heating could also be ad- 
dressed analytically. The spherical collapse model, in which 
the collapsing object is divided into a series of concentric 
shells, seems better suited than the ellipsoidal collapse. A 
thorough discussion of tidal heating is beyond the scope of 
the present paper. Note, however, that the torque imparted 
by the (external) mass distribution on a thin spherical shell 
is r(x) oc fdttS(x) A V<j>{x) (e.g. Ryden & Gunn 1987), 
where 5(x) — 8(x) — 5(x) is the deviation from the spheri- 
cally symmetric distribution S(x). This is 5 which pulls the 
infalling matter out of its purely radial motion. In the lin- 
ear regime, S is a Gaussian density field statistically inde- 
pendent of the spherical average S(x). Hence, tidal heating 
does not induce any environmental dependence at first or- 
der. Consequently, modelling the growth of nonlinearities 
in the surrounding mass distribution will be crucial to as- 
certain analytically the importance of tidal heating in the 
environmental dependence of halo collapse. 



6 CONCLUSION 

The ellipsoidal collapse model is an extension of the spheri- 
cal dynamics that takes into account the anisotropic collapse 
of triaxial perturbations. This non-spherical dynamics pro- 
vides a substantially better description of halo statistics such 
as mass function and large scale bias. It is, however, unclear 
whether the ellipsoidal collapse can induce environmental 
effects similar to that seen in N-body simulations. 

In this paper, we have attempted to address this issue, 
paying special attention to both the statistical and dynam- 
ical origin of the environmental dependence. In a first part, 
we have explored the statistical correlation that arises in 
(Gaussian) initial conditions between the local properties of 
the shear and the configuration of the large scale environ- 
ment. To this purpose, a number of joint statistics for the 
shear tensor have been derived, thereby extending the pre- 
vious analysis of Doroshkevich (1970). In a second part, we 
have examined the dynamical aspect of the environmental 
dependence using a simplified model that takes into account 
the interaction between a collapsing, ellipsoidal perturbation 
and its large scale environment. Relaxing the assumption 
of sphericity (at the heart of the spherical collapse) intro- 
duces a dependence of collapse redshift on environment. The 
tidal force exerted by the surrounding mass distribution al- 
ters the collapse of the major and minor axes, and causes 
haloes embedded in large overdensities to virialize earliest. 
We have found that the environment density is a key pa- 
rameter in determining the virialization redshift, the large 
scale asphericity contributing mostly to increase the scatter 
in collapse density. An effective barrier whose shape depends 
on the large scale density provides a good description of this 
environmental effect. Such an interpretation has the advan- 
tage that the EPS formalism can be applied to estimate the 
environmental dependence of halo properties like formation 
redshift and large scale bias. 

We have shown that, using this moving barrier ap- 
proach, a correlation between formation redshift, large scale 
bias and environment density naturally arises. The magni- 
tude of the effect is similar, albeit smaller, to that seen in N- 
body simulations. It is large for low mass haloes M <C M*, 
and fades as we go to high masses M > M+ as a result 
of a genuine statistical effect, namely, the decrease in aver- 
age asymmetry and stochasticity with increasing halo mass. 
Haloes that formed at high redshift are found to be more 
clustered than those that assembled recently. This is pre- 
cisely the behaviour reported by Sheth & Tormen (2004); 
Gao, Springel & White (2005). However, haloes in denser 
regions are predicted to assemble later. This result is incon- 
sistent with the trend measured in overdense environments 
S J>0. It calls into question the role of the ellipsoidal collapse 
in shaping the halo mass function. Nevertheless, several lines 
of recent evidence indicate that, in relatively underdense 
regions, the average formation redshift increases with de- 
creasing environment density (e.g. Harker et al. 2006), while 
haloes may be more strongly clustered than the mean (Ab- 
bas & Sheth 2007). Our model predicts and effect in the right 
sense. This suggests that the ellipsoidal collapse model may 
be applicable in underdense regions, where tidal interactions 
are weak or absent. Conditional halo mass functions n(M\S) 
could provide another testable prediction since, in the mov- 
ing barrier interpretation adopted here, it is expected that 
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n(M\S) in underdense regions should be (slightly) biased 
towards high mass haloes as compared to the prediction of 
Sheth & Tormen (2002). 

Recently, Sandvik et al. (2007) have discussed a multi- 
dimensional extension of the EPS formalism that takes into 
account ellipsoidal collapse (Chiueh & Lee 2001). They 
find a very weak correlation between halo assembly his- 
tory and environment, presumably because their implemen- 
tation includes only the statistical aspect of environmen- 
tal effects. In our scenario, the dynamical interaction be- 
tween the external mass distribution and the collapsing 
halo plays a crucial role in the environmental dependence 
of halo properties. The statistical correlations contribute 
mostly to increase the strength of the effect with decreas- 
ing halo mass. Keselman & Nusser (2007) have shown that 
the age-dependence of halo bias persists in a simplified de- 
scription of gravitational dynamics, suggesting thereby that 
(at least some of) the environmental dependence arises in 
the early stages of the collapse. They find that halo for- 
mation redshift strongly correlates with the dimensionless 
parameter -q oc (1 + &e\ + 2pi)~ 1 ^ 2 , young haloes forming 
from fluctuations with higher r\ than old ones. They argue 
that this follows from the dependence of first axis collapse 
on configuration shape, namely, a planar perturbation col- 
lapses faster than a spherical perturbation with the same 
initial density (Bertschinger & Jain 1994) . Their results are, 
however, difficult to reconcile with the above interpretation 
since it implies that, for a fixed halo mass, haloes forming in 
regions of larger r\ are older (because they grow less rapidly). 
This discrepancy can be alleviated if one assumes that the 
formation of virialized haloes corresponds to the collapse of 
the third axis. In this case, the collapse is always faster for 
spherical perturbations (Audit, Teyssier & Alimi 1997) so 
that young haloes tend to form in regions of higher rj, i.e. 
more spherical or denser environments, in agreement with 
our findings. 

A serious shortcoming of our model is the neglect of 
anisotropies beyond the quadrupole term in the external 
density field. Furthermore, apart from a global rotation, it 
ignores the non-radial degrees of freedom in the collapsing 
protohalo. It would be of great interest to ascertain whether 
non-radial motions created by a clumpy, growing large scale 
mass distribution can induce an effect similar to that seen 
in simulations. A priori, tidal heating should be more effi- 
cient in high density environments. It may plausibly reverse 
the trend found in this paper, namely, increase the critical 
collapse density in large scale overdensities. If this happens 
to be true, a moving barrier approach would naturally pre- 
dict haloes in dense regions to assemble relatively early and 
to be more strongly clustered than the mean. Alternatively, 
Wang et al. (2007) have suggested that tidal heating causes 
haloes to appear less massive than expected from their initial 
density field. This may also provide an explanation for the 
relatively large/low bias of old/young haloes. Clearly, large 
N-body simulations are needed to understand the complex 
relations reported by, e.g. Harker et al. (2006), Wechsler et 
al. (2006), Gao & White (2007). This caveat notwithstand- 
ing, analytic models can provide an elegant route to captur- 
ing the essential features of the environmental dependence. 

In standard galaxy formation models, the correlation 
between the haloes and their large scale environment intro- 
duces a correlation between the properties of galaxies and 



the regions they occupy (e.g. Navarro, Abadi & Steinmetz 
2004; Abbas & Sheth 2005, 2006; Berlind et al. 2005). The 
observational results of Skibba et al. (2006), Blanton et al. 
(2006) and Tinker et al. (2007) support this prediction and 
leave little room for a galaxy assembly bias (see, however, 
Croton et al. 2007). Nevertheless, the influence of the en- 
vironmental dependence of halo properties on the galaxy 
population remains unclear. It would be valuable to assess 
whether environmental effects produced by the anisotropic 
collapse of haloes can leave a detectable signature in the 
properties of field galaxies. 
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APPENDIX A: TWO-POINT CORRELATIONS 
OF THE SHEAR TENSOR 

We consider the two-point correlation functions of an arbi- 
trary symmetric tensor field Tij(x). Statistical isotropy and 
symmetry imply that, in position space, these correlations 
must be of the form 

(Ty(x)T !m (x+ r)) = ^li^fifjhfm 

+ ^2(r) {hnSjrn. + fif m 5jl + fjflSim + fjf m Su) 

+ * 5 (r) (SuSjm + S im 5ji) , (Al) 

where fi = Ti/\r\ and the functions \&i(r) depend on r = 
\r\ only. This is the most general ansatz for a symmetric, 
isotropic correlation tensor. In the case of a scalar (spin-0) 
tensor such as the shear fjj defined in eq. fl}, $2 = ^3 and 
^4 = Explicit expression for the functions &i(r) can be 
obtained from the Fourier transform of {£ij(x)£i rn (x + r)}, 



(Cij( fc )Clm(fc)) = Pg(k) kikjklkr, 



(A2) 



where ki = ki/k and Ps(k) is the power spectrum of the 
density field S (x). The following integrals are useful to the 
calculation of 9i(r), 



-jo(kr) - -j 2 (fcr) 



1 / 

^ I i 2 ikrii 

- I d/i/i e 

1 r +1 

- J dfj.fi e 



Tjo(kr) 



3 

j32(kr) + ^Mkr) , (A3) 



where je (x) are spherical Bessel functions of the first kind. 
With these informations, the functions 'J; may be conve- 
niently expressed as 



*i(r) = 
*s(r) = 
*s(r) = 



dlnfc A 2 S (k) ji(kr) 



(A4) 



dlnfeA|(fe) 



dlnfcAj(fc) 
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Figure Al. The correlation functions \Pi(r), ^(r) and ^s(r) 
as a function of comoving separation r for the ACDM cosmology 
considered in this paper. They are normalised by £01 ( e 1- HOD - 
The shear tensor at position x and x + r is smoothed on scale 
Ro = 10 and R% = 2 ft _1 Mpc, respectively. Note that is 
strongly suppressed on scale r <,Ro- 



into the form given in, e.g., Lee & Pen (2001), Crittenden 
et al. (2001), and Catelan & Porciani (2001), which involves 
J n = nr~ n f Q ds £(s)s n ~ 1 (this is best seen by integrating J n 
by part). 



APPENDIX B: AN INTEGRAL OVER THE 
SO(3) MANIFOLD 

In this Appendix, we present the calculation of the integral 
eq. (|19[) over the special orthogonal group SO(3), where a 
and A are real symmetric matrices. 

Let R = R fl (</j) be the rotation by an angle ip about the 
axis h. Since the Haar measure dR is invariant under both 
left and right transformations, we can redefine R so as to 
diagonalise A and a. We will therefore assume that the ma- 
trices A and a are diagonal, A = diag(A^) and a — diag(Qi). 
It is worth noting that a and A are not orthogonal matrices, 
but belong to the (real) linear group GL(3,R). Consider now 
the rotation of 7r around the first axis, o\ = diag(l, —1, — 1), 
and its distinct permutations 02 = diag(— 1, 1, — 1) and 
0-3 = diag(— 1, — 1, 1). Since the compositions cr^R char- 
acterise equivalent principal axis frames, the integral (|19|) 
should in principle be performed over the quotient space 
SO(3)/{<7i, 1T2, as}. However, the OiS are diagonal and we 
find 



tr [((TjR) A (cr l R) T a] = tr (RAR T a) 



(Bl) 



In the limit r — > 0, both and ^3 vanish, but ^5 = £oi/15 
(Figure [Al|) . Notice that the functions can also be cast 



Therefore, one can also perform the integration over the 
whole orthogonal group SO (3) and multiply the final result 
by 1/4. 

Integrals of the form 

F(a,X) = fdfie^f^^ 1 ") (B2) 
Jg 

where G (fi GG) is a compact Lie group and a, A belong 
to the linear group GL, appear in Quantum Field Theory, 
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chemistry as well as in harmonic analysis (e.g. Hua 1963; 
Itzykson & Zuber 1980; Wei & Eichinger 1989). When G is 
the unitary group U(N), group theoretical techniques such 
as character expansion can be used to obtain a closed form 
determinantal evaluation of (IB2|I (e.g. Harish Chandra 1958; 
Balantekin 2000). Unfortunately, these methods cannot be 
easily applied when G is the orthogonal group SO(N), es- 
sentially because the representations of SO(N) and GL(N,R) 
are very different. Instead, one generally relies on a specific 
parametrisation of the rotation matrices and express the re- 
sult in terms of a series of orthogonal functions. 

Here we follow the approach outlined in Wei & Eichin- 
ger (1990) and parametrise the special rotation matrices R 
in terms of the Euler angles < (p, ip < 2n, < $ < n, 

— s^c v — c$s v c^ —s^s v + c$c v c$ c,ps^ , (B3) 

where = costp, s v — sirup etc. Notice that this represen- 
tation becomes singular when t? = or n (such singularities 
are expected owing to the topology of SO (3)). The origi- 
nal integral (|19[) over the SO(3) manifold can be written 
as a triple integral with the (normalised) invariant measure 
dR= l/(87r 2 )sinj?dipdi9di/), 

F(a,\) = —r d<p dtfsintf / dtp e^^^ . (B4) 
87r Jo Jo Jo 

Since the integrand is invariant under the transformations 
(a, A) -» (fca,|A) and (a, A) - ( Q + fcl, A - l), 
where k is a real number and I is the 3x3 identity ma- 
trix, F(a, A) can be expressed as a function of four variables 
instead of the original six matrix eigenvalues. This can also 
seen by rewriting the trace as (Wei & Eichinger 1990) 



tr (RAR T a) = - (tra trA - tra trA) + tr (RAR 7 



(B5) 



where, for instance, a — diag («i3, a 2 3, 0), A = 
diag (Ai 3 , A23, 0), onj = on — otj and A,j = A; — Xj. We 
enforce the ordering ai > a 2 > «3 and Ai > A2 > A3, 
so that ai3,a 2 3 > and Ai 3 ,A 2 3 > 0. The integral de- 
pends now on four distinct combinations of the as and 
As: e+ = (l/3)tratrA, e_ = (l/3)tratrA, e a = cxiz/tia, 
ex = Ai 2 /trA. With this parametrisation, —00 < e + < 00, 
e_ > and < e a , ex < 1. 

We can expand the trace in terms of the Wigner D- 
functions T> m , I being the index of the representation 
(see, e.g., Sakurai 1985). These 3D harmonics generate ir- 
reducible representations of the three-dimensional rotation 
group and, therefore, form a complete orthogonal set of 
functions defined on SO (3) itself. Unsurprisingly, only the 
quadrupole (I = 2) rotation matrices appear in the trace 
decomposition, 

"2 



tr RAR a 



3 

2 £a 



V 

3 0M 



+ \^e x (v\ 



■v: 



>+K-*) 



(B6) 



The explicit form of these I = 2 harmonics is given in Ta- 
ble [Bl] Note that tr(RAR T a) depends on the three "shape" 
parameters e_, e a and ex solely, because points on SO(3) 



truly have only three degrees of freedom. The integral over 
the variable ip can be performed using the following identity 
(see §3.937 of Gradshteyn & Ryzbik 2000), 



2tt 



2n Jo V 



(B7) 



where lo is a modified Bessel function of the first kind. The 
result can be written 



F{a,X) = e l3e +W(l3e-,e a ,e x ) 



(B8) 



where the function W (/3e_ , e a , e\) is a double integral which 
can be arranged such that 



W{Pe_,e a ,e x ) 

-3 e _ I 1 



X/l 



2tt 

3/3e_£A 



dr 



dip exp 



3/3 e- 



g(r,ip,e a ) 



\J h(r, <p, e a ) j , 



(B9) 



where r = — cost). The functions g(r,ip,e a ) and h(r,ip,e a ) 
are defined as 

g(r, <p,e a ) = 1 + r 2 + e a (l - r 2 ) cos(2ip) 

h(r,<p,e a ) = g 2 - 4 (l - e 2 ) r 2 . (B10) 

They are periodic of period tt in the argument (p. Further- 
more, on the domain defined by < r < 1 and < tp < 2n, 
the function g is bounded by 1 — e a < g(r, ip,e a ) < 1 + £<*• 
Consequently, the double integral that appears in eq. (|B9|I 
is always larger than or equal to 1. The equality holds only 
when a and/or A is proportional to the identity matrix. 

Expanding the integrand of ()B9|I about /3e_ = gives, 
upon integration, the following series 



W {pe- 



ea, ex) 
I3 2 e 2 _ 



(Bll) 



1 + 



40 



(3 3 e 3 _ 

px(e a )pi(ex) + -— 7 -p2{e a )p2{e x ) 



+ 



4480 Pl(£a)Pl(£A)+ 73920 



840 

Pi {e a )p2 (e a )pi (ex )p 2 (ea) 



where pi(x) — (Zx 2 + 1) and P2(x) = (9x 2 — 1). We have 
found that this truncated, fifth-order expansion is accurate 
to within 2 per cent in the range < /3e_ <T-5 and can be 
used efficiently in the computation of the group integral (1191) . 
Conditional probability distributions for the relative orien- 
tation can be derived in a straightforward way from the 
preceding results. 

Noteworthy is the special case G=SO(2), in which the 
integral (|B2|) can be written in closed form, 



F(a,X) = e /3E +/ (~ai 2 Ai 2 ) , (B12) 
where the parameter e+ is now e+ — (1/2) tra trA. 



APPENDIX C: POTENTIAL OF A 
HOMOGENEOUS ELLIPSOID 

For a homogeneous ellipsoid defined such that p(r) = Sp if 
r T (AA T ) -1 ) r< 1, and zero otherwise, the potential at any 
exterior point r = (n, r 2 , 7-3) is (Kellog 1929; Chandrasekhar 
1969) 
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Table Bl. Quadrupole Wigner D-functions T> 2 (i/?,i9,-0) (in ZXZ representation. Harmonics with mi, 7712 = ±1 are not shown). 







m 2 = —2 


m 2 = 


m 2 = 2 


fill 


= -2 


i(l + cosi?) 2 e 2 ^+ 2 "'' 


-^/fsin^e 2 ^ 


i (1 -cosi?) 2 e 2 ^- 2 ^ 


mi 


= 


-n/fsin^e 2 ^ 


\ (3cos 2 i? - l) 


-^/fsi^iSe" 2 ^ 


mi 


= 2 


i (1 - cost?) 2 e - 2i f+ 2i * 


-^/|sin 2 i?e- 2 ^ 


1 (1 + cosi?) 2 e- 2 ^- 2i1 '' 



$ e (r) = 7rG5p 



(CI) 



while the potential inside and on the boundary of the ellip- 
soid is 



<!>i(r) = nGSp 



(C2) 



AA is a positive definite matrix whose eigenvalues are the 
square of the principal axis lengths A^ > A2 > A- A > 0. The 
variable p is defined as the algebraically largest root of the 
following cubic equation, 

222 
x y z _ 

+ + — = 1 



A{ + p Ai+p Aj + p 
The functions C(p) and bi(p) are given by 
du 



(C3) 



C(p) = AiA 2 A 3 
bi(p) = A!A 2 A 3 



A(«) 



L4f + u)A( W ) 



(C4) 



where A(it) = n| =1 LA 2 , + w) '". These functions can be 
expressed in terms of the Legendre's incomplete elliptic in- 
tegrals of first and second kind (Binney & Tremaine 1987). 
The gravitational potential energy of this self-gravitating el- 
lipsoid is 



Spot = — yjj GM 2 



du 



(C5) 



and is proportional to the potential at x = 0, C(0). Notice 
also that Poisson's equation implies ^\ 6j(0) = 2. Eq. ()C5|) 
corrects a typographical error in eq. (24) of Eisenstein & 
Loeb (1995). 
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